1*1 


Defence  Research  and  Recherche  et  developpement 
Development  Canada  pour  la  defense  Canada 


Copy  No. 


Modelling  and  Simulating 
Unsteady  Six  Degrees-of-Freedom 
Submarine  Rising  Maneuvers 


George  D.  Watt 


Defence  R&D  Canada  -  Atlantic 

Technical  Report 
DRDC  Atlantic  TR  2007-008 
February  2007 


Canada 


This  page  intentionally  left  blank. 


Modelling  and  Simulating 
Unsteady  Six  Degrees-of-Freedom 
Submarine  Rising  Maneuvers 


George  D.  Watt 


Defence  R  &  D  Canada  -  Atlantic 

Technical  Report 

DRDC  Atlantic  TR  2007-008 

February  2007 


Author 


Original  signed  by  George  D.  Watt 

George  D.  Watt 


Approved  by 

Original  signed  by  Neil  Pegg 

Neil  Pegg 

Head,  Warship  Performance 


Approved  for  release  by 


Original  signed  by  J.L.  Kennedy 


K.  Foster 

Chair,  Document  Review  Panel 


©  Defence  Research  and  Development  Canada,  2007 
©  Recherche  et  developpement  pour  la  defense  Canada,  2007 


Abstract 


DRDC  Atlantic  is  collaborating  with  ANSYS  Canada  and  the  University  of  New  Brunswick  to 
develop  an  unsteady,  six  degrees-of-freedom,  Reynolds  Averaged  Navier-Stokes  (RANS)  sub¬ 
marine  maneuvering  simulation  capability.  Initially,  this  is  being  used  to  evaluate  emergency 
rising  maneuvers.  During  these  maneuvers,  high  negative  angles  of  attack  can  occur  which 
result  in  a  roll  instability  not  previously  predicted  by  quasi-steady  modelling.  The  objective  of 
the  RANS  simulation  is  to  reproduce  the  instability  and  investigate  mitigation  strategies. 

Models  for  weight  and  buoyancy  when  blowing,  high  incidence  propulsion,  and  appendage 
and  propulsion  activation  are  presented  and  tested.  A  high  incidence,  quasi-steady,  coeffi¬ 
cient  based  hydrodynamic  model  used  in  previous  stability  analyses  is  also  presented.  These 
models  are  used  for  evaluating  stability,  testing  the  system  models,  and  investigating  different 
maneuvering  scenarios  in  preparation  for  carrying  out  the  computationally  intensive  RANS 
simulations.  These  preliminary  investigations  suggest  the  initial  roll  angle  prior  to  blowing 
ballast,  coupled  with  the  roll  instability  and  low  pitch  angles,  plays  an  important  role  in  the 
emergence  roll  angle. 


Resume 

RDDC  Atlantique,  en  collaboration  avec  ANSYS  Canada  et  l’Universite  du  Nouveau-Brunswick 
(UNB),  developpe  une  capacite  de  simulation  de  manoeuvres  instables  d’un  sous-marin  a  six 
degres  de  liberte  par  l’application  d’equations  de  Navier-Stockes  a  moyenne  de  Reynolds  (RANS). 
Initialement,  cette  capacite  est  utilisee  pour  evaluer  les  manoeuvres  de  remontee  d’urgence.  Au 
cours  de  ces  manoeuvres,  des  angles  d’attaque  negatifs  eleves  peuvent  etre  obtenus,  entrainant 
une  instability  du  roulis  qui  n’etait  pas  prevue  anterieurement  par  la  modelisation  quasi-stable. 
L’objectif  de  la  simulation  RANS  est  de  reproduire  l’instabilite  et  d’en  etudier  les  possibility 
d’ attenuation. 

Les  modeles  de  poids  et  de  flottabilite  lors  du  delestage,  la  propulsion  a  haute  incidence  et 
1’ activation  de  l’appendice  et  de  la  propulsion  sont  presentes  et  verifies.  On  presente  egalement 
un  modele  hydrodynamique  de  haute  incidence,  quasi-stable  et  base  sur  les  coefficients  qui 
est  utilise  dans  des  analyses  anterieures  de  stabilite.  Ces  modeles  sont  utilises  pour  evaluer 
la  stabilite,  verifier  les  modeles  de  systemes  et  etudier  les  differents  scenarios  de  manoeuvre  en 
preparation  a  Pexecution  de  simulations  exigeant  un  grand  nombre  de  calculs  RANS.  Les  etudes 
preliminaires  revelent  que  l’angle  de  roulis  initial  avant  le  delestage,  jumele  a  l’instabilite  du 
roulis  et  aux  faibles  angles  de  tangage,  joue  un  role  important  dans  l’angle  du  roulis  d’emersion. 
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Executive  Summary 


Introduction 

DRDC  Atlantic  is  collaborating  with  ANSYS  Canada  and  the  University  of  New  Brunswick 
(UNB)  to  develop  an  unsteady,  six  degrees-of-freedom  (DOF),  Reynolds  Averaged  Navier-Stokes 
(RANS)  simulation  of  a  maneuvering  submarine.  This  capability  is  initially  being  used  to 
evaluate  emergency  rising  maneuvers.  In  a  typical  rising  scenario,  ballast  is  blown  at  depth, 
the  nose  of  the  boat  is  pitched  up,  and  speed  is  increased  to  minimize  large  negative  angles 
of  attack.  Despite  these  precautions,  full  scale  trials  show  that  an  underwater  roll  instability 
still  can  occur  and  can  result  in  excessive  roll  when  the  boat  surfaces.  Quasi-steady  modelling 
has  not  previously  predicted  satisfactorily  the  roll  angles  resulting  from  the  instability.  The 
objective  of  the  collaboration  is  to  develop  the  unsteady  RANS  simulation  capability,  use  it  to 
reproduce  the  instability  and,  if  successful,  investigate  mitigation  strategies. 

This  report  presents  the  6  DOF  equations  of  motion,  a  coefficient  based  quasi-steady 
hydrodynamic  model  (a  stand-in  for  the  RANS  model),  and  several  system  models  required  to 
support  the  simulations.  It  uses  these  models  to  investigate  various  rising  scenarios,  showing 
how  blowing,  the  sternplanes,  propulsion,  control,  and  the  initial  heel  angle  of  the  boat  effect 
the  rising  maneuver. 

Significance 

The  unsteady  RANS  simulation  capability  is  a  step  towards  improving  our  ability  to  predict 
and  understand  operational  limitations.  It  can  be  used  in  many  maneuvering  scenarios.  The 
system  models  presented  (blowing,  control,  and  propulsion  models)  are  used  in  both  the  RANS 
and  coefficient  based  simulations.  The  coefficient  based  simulations  provide  an  efficient  test 
bed  for  preliminary  investigations  that  will  minimize  the  number  of  computationally  intensive 
RANS  simulations  required  later. 

Principal  Results 

The  models  are  described,  successfully  implemented,  and  preliminary  rising  simulations  carried 
out  using  the  coefficient  based  hydrodynamic  model.  This  exercises  the  models  and  establishes 
suitable  scenarios  for  the  RANS  simulations.  Unexpectedly,  this  preliminary  work  has  shown 
that  initial  heel,  together  with  the  roll  instability  and  low  pitch  angles,  significantly  increases 
roll  angles  while  rising. 

Further  Work 

ANSYS  Canada  has  completed  development  of  the  unsteady  6  DOF  RANS  simulation  capa¬ 
bility.  The  capability  is  being  transferred  to  UNB  where  rising  simulations  will  be  carried  out 
and  compared  with  the  coefficient  based  simulations  presented  herein.  This  will  provide  an 
opportunity  to  validate  the  current  results  and  further  investigate  the  role  of  the  initial  heel 
angle.  This  work  should  be  supplemented  by  full  scale  trials  to  see  what  heel  angles  occur  in 
practice  and  to  see  if  mitigation  strategies  suggested  by  the  simulations  are  realistic. 


G.D.  Watt,  2007.  Modelling  and  Simulating  Unsteady  Six  Degrees-of-Freedom  Submarine 
Rising  Maneuvers.  DRDC  Atlantic  TR  2007-008.  Defence  R&D  Canada  -  Atlantic. 
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Sommaire 

Introduction 

RDDC  Atlantique,  en  collaboration  avec  ANSYS  Canada  et  l’Universite  du  Nouveau-Brunswick, 
developpe  une  capacite  de  simulation  de  manoeuvres  instables  d’un  sous-marin  a  six  degres 
de  liberte  par  l’application  d’equations  de  Navier-Stockes  a  moyenne  de  Reynolds  (RANS). 
Initialement,  cette  capacite  est  utilisee  pour  evaluer  les  manoeuvres  de  remontee  d’urgence.  Dans 
un  scenario  typique  de  remontee,  le  ballast  est  deleste  en  profondeur,  le  devant  du  sous-marin 
est  deplace  vers  le  haut  et  la  vitesse  est  accrue  pour  reduire  les  angles  d’attaque  negatifs  eleves. 
Malgre  ces  precautions,  les  essais  complets  demontrent  que  le  roulis  sous-marin  peut  encore 
etre  instable  et  atteindre  un  niveau  excessif  lors  de  l’emersion  du  sous-marin.  La  modelisation 
quasi-stable  ne  permet  pas  encore  de  predire  de  fagon  satisfaisante  Tangle  de  roulis  cause  par 
Pinstabilite.  L’objectif  de  la  collaboration  est  de  mettre  sur  pied  la  capacite  de  simulation 
RANS  instable,  de  l’utiliser  pour  reproduire  Pinstabilite  et,  le  cas  echeant,  d’en  etudier  les 
possibility  d’ attenuation. 

Le  present  rapport  presente  les  equations  de  mouvement  a  six  degres  de  liberte,  un  modele 
hydrodynamique  quasi-stable  base  sur  les  coefficients  (remplagant  le  modele  RANS),  et  plusieurs 
autres  modeles  requis  pour  appuyer  les  simulations.  Le  rapport  utilise  ces  modeles  pour  exposer 
differents  scenarios  de  remontee  afin  de  demontrer  comment  le  ballast,  les  tableaux  arriere,  la 
propulsion,  le  controle  et  Tangle  de  gite  initial  du  sous-marin  ont  un  impact  sur  la  manoeuvre 
de  remontee. 

Portee 

La  capacite  de  simulation  RANS  quasi-stable  est  une  etape  vers  P amelioration  de  notre  habilete 
a  predire  et  comprendre  les  limites  operationnelles.  Cette  capacite  peut  etre  utilisee  dans  de 
nombreux  scenarios  de  manoeuvres.  Les  modeles  de  systeme  presentes  (modeles  de  delestages, 
de  controle  et  de  propulsion)  sont  utilises  dans  les  simulations  basees  sur  RANS  et  les  coeffi¬ 
cients.  Les  simulations  basees  sur  les  coefficients  fournissent  un  banc  d’essai  eventuel  pour  des 
recherches  preliminaires  qui  reduiront  les  simulations  a  grand  nombre  de  calculs  RANS,  qui 
seront  necessaires  plus  tard. 

Result  at  s 

Les  modeles  sont  decrits  et  appliques  avec  succes,  et  des  simulations  preliminaires  de  remontee 
sont  effectuees  a  l’aide  d’un  modele  hydrodynamique  base  sur  les  coefficients.  Les  modeles  sont 
ainsi  mis  a  l’epreuve  et  servent  a  etablir  des  scenarios  adaptes  aux  simulations  RANS.  Contre 
toute  attente,  ce  travail  preliminaire  a  demontre  que  le  gite  initial,  jumele  a  Pinstabilite  du 
roulis  et  aux  faibles  angles  de  tangage,  augmente  de  fagon  significative  les  angles  du  roulis  lors 
de  la  remontee  du  sous-marin. 

Recherches  futures 

ANSYS  Canada  a  termine  le  developpement  de  la  capacite  de  simulation  RANS  instable  a  six 
degres  de  liberte.  La  capacite  est  transferee  a  l’UNB  ou  des  simulations  de  remontee  seront 
effectuees  et  comparees  avec  les  simulations  basees  sur  les  coefficients  qui  sont  presentees  ici.  II 
sera  done  possible  de  valider  les  resultats  actuels  et  de  pousser  l’etude  sur  le  role  de  Tangle  de  gite 
initial.  Ce  travail  devrait  etre  enrichi  par  des  essais  complets  visant  a  determiner  les  angles  de 
gite  reels  et  a  verifier  si  les  possibility  d’attenuation  suggerees  par  les  simulations  sont  realistes. 

G.D.  Watt,  2007.  Modelisation  et  simulation  de  manoeuvres  instables  de  remontee  d’un  sous- 
marin  a  six  degres  de  liberte.  RDDC  Atlantique  TR  2007-008.  R&D  pour  la  defense  Canada 
-  Atlantique. 
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Nomenclature 


CB,  CG 

d 

D 

B  =  pV  g 
DG  =  zG  —  zB 

g 

I 

J,  J}, 

K,  M,  N 
K',M',N' 

KP 

Kt-Kq 

l 

m,  m0 

MBT 

n 

N 

Pa 

p,q,r 

Q 

R 

R  =  Ul/v 
t,  te 
t. 

T 

T 

U ,  V,  w 

wT,wT0 

U  =  \/u2  +  v2  +  w2 

Us 

V 


Centers  of  buoyancy  and  gravity. 

Maximum  hull  diameter. 

Propeller  diameter. 

Buoyancy. 

Height  of  the  CB  above  the  CG. 

Gravitational  constant. 

Moments  of  inertia  in  body  axes. 

Propeller  advance  ratio  and  behind-the-boat  advance  ratio. 
Body  axis  moments. 

Body  axis  moments  nondimensionalized  by  pU2l3 /2. 

Effect  of  propeller  torque  on  rolling  moment  K . 

Propeller  thrust  and  torque  coefficients. 

Overall  length  of  the  hull. 

Mass,  initial  mass  within  V . 

Main  Ballast  Tank;  several  at  different  axial  locations. 

Propeller  revolutions  per  second,  rps. 

Number  of  MBTs. 

Atmosperic  pressure. 

Body  axis  angular  velocities. 

Propeller  torque. 

Gas  constant  for  air. 

Reynolds  number. 

Time,  time  of  emergence. 

Propeller  thrust  deduction  (in  §4). 

Thrust  (in  §4). 

Temperature  (in  §5). 

Body  axis  velocities. 

Taylor  wake  fraction,  generally  and  at  zero  incidence. 

Overall  speed  of  vehicle. 

Roll  stability  index,  m/s. 

Volume  of  the  external  hydrodynamic  envelope,  including  main 
ballast  tanks. 
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Va,Vai 

VA 

VT,VTi 
W  =  mg 
x,y,z 
*To>  Vo  >  zo 

xBiVBi zB 
XG1  Vg-i  zG 
XT  i  1  ZT  i 

xn,z» 

X,  Y,  Z 

X',Y',Z' 

XP 


Volume  of  air  in  all  MBTs,  and  in  tank  i. 

Speed  of  advance  of  propeller. 

Volume  of  all  MBTs,  and  of  tank  i. 

Weight  within  V. 

Body  fixed  axes. 

Inertial  (earth-fixed)  axes. 

Coordinates  of  CB  (centroid  of  V)  in  body  axes. 

Coordinates  of  CG  (center  of  mass  m )  in  body  axes. 
Coordinates  of  the  centroid  of  MBT  i ;  yTi  =  0 . 

Coordinates  of  the  centroid  of  the  blown  mass  fraction;  =  0 . 
Body  axis  forces. 

Body  axis  forces  nondimensionalized  by  pU2£2 /2. 

Effect  of  propeller  thrust  on  axial  force  X. 

Depth  below  the  ocean  surface  of  the  water  level  in  MBT  i. 


a  =  tan  1(w/u) 

/ 3  =  tan  ~1{—v/u) 

^  b 

0  =  tan_1(\/u2  +  w2/u ) 

v 

P 

$  =  tan_1(— vj— w) 


Angle  of  attack. 

Angle  of  drift. 

Rudder,  stern,  and  foreplane  deflections;  direction  is  found  from 
the  right  hand  rule  using  body  axes. 

Flow  incidence,  always  positive. 

Blown  mass  fraction,  overall  and  for  tank  i. 

Kinematic  viscosity  of  water. 

Density  of  sea  water. 

Flow  orientation. 

Yaw,  pitch,  and  roll  Euler  angles  giving  body  axes  orientation 
relative  to  inertial  axes. 

Initial  heel,  the  roll  angle  prior  to  blowing  ballast. 


Subscripts 

a  The  blown  air  in  the  MBTs. 

c  A  command. 

o  The  condition  immediately  prior  to  blowing  the  MBTs. 

t  The  trim  component. 

T  Main  ballast  tank. 

U  Value  at  which  the  vehicle  becomes  unstable  (when  U$  =  0). 

A  dot  over  a  symbol  indicates  differentiation  with  respect  to  time;  eg,  u  =  du/dt.  The  following 
relations  are  useful: 

u  =  U  cos  0  v  =  —U  sin  0  sin  <!> 

w  =  —U  sin  0  cos  $  \/ v2  +  w2  =  U  sin  © 


vii 


DRDC  Atlantic  TR  2007-008 


Acknowledgment 

The  author  very  much  appreciated  feedback  from  Mr.  Mark  Bettle,  graduate  student,  University 
of  New  Brunswick,  which  enabled  important  corrections  to  be  made  to  this  report. 


DRDC  Atlantic  TR  2007-008 


1  Introduction 


6  DOF  RANS  Simulation 

DRDC  Atlantic  is  collaborating  with  ANSYS  Canada  and  the  University  of  New  Brunswick 
(UNB)  to  develop  an  unsteady,  six  degrees-of-freedom  (DOF),  Reynolds  Averaged  Navier-Stokes 
(RANS)  submarine  maneuvering  simulation  capability.  This  capability  is  needed  to  validate 
fast,  coefficient  based  simulations  used  to  investigate  maneuvering  limitations  and  establish  safe 
operating  envelopes  for  underwater  vehicles.  Several  countries  (eg,  the  US,  UK,  France)  use 
a  free  swimming  scale  model  for  such  validations  which,  according  to  the  US  [l],  “is  currently 
the  best  predictor  of  full  scale  submarine  maneuvering  performance.”  These  facilities  cost  tens 
of  millions  of  dollars  to  develop,  maintain,  and  use.  A  computational  fluid  dynamics  (CFD) 
validation  is  much  cheaper  and  capable  of  providing  better  detail.  The  disadvantage  to  a  CFD 
capability  is  that  its  predictions  are  not  as  reliable  as  experiments.  But  CFD  technology  is 
evolving  quickly  so  it  is  worth  developing  and  evaluating  this  capability.  By  collaborating 
with  a  successful  commercial  CFD  vendor,  there  is  potential  for  commercializing  the  capability 
thereby  minimizing  ongoing  maintenance  and  development  costs. 

ANSYS  Canada  has  developed  a  basic  capability  for  use  with  its  commercial  RANS  code 
CFX  [2].  It  requires  that  the  flow  field  be  discretized  with  a  rigid,  body  fixed  mesh  extending 
from  the  surface  of  the  vehicle  out  to  the  far  field.  The  mesh  and  boat  move  together  controlled 
by  the  same  6  DOF  solid  body  equations  of  motion  used  by  the  DRDC  Submarine  Simulation 
Program  (DSSP)  [3].  Unlike  DSSP,  which  inputs  quasi-steady  hydrodynamic  information  and 
solves  the  equations  of  motion,  CFX  solves  the  unsteady  RANS  equations  for  the  fluid  flow  about 
the  submarine,  calling  the  solid  body  equations  of  motion  for  information  on  how  to  change 
the  flow  boundary  conditions  at  each  time  step.  CFX  passes  the  unsteady  hydrodynamic  forces 
to  the  solid  body  equations  which  account  for  inertia,  buoyancy,  propulsion,  and  control  forces 
using  the  models  presented  below.  The  CFX  simulation  is  at  best  second  order  accurate  in 
space  and  time  while  DSSP  is  fourth  order  accurate  in  time.  However,  the  CFX  time  steps  are 
much  smaller  than  those  used  by  DSSP  since  they  are  determined  by  the  complex,  unsteady 
hydrodynamic  flow  being  modelled. 

Currently,  the  CFX  simulation  can  model  any  maneuver  in  which  the  boat  is  deeply  sub¬ 
merged  and  isolated  from  any  other  vehicle  or  boundary.  In  developing  the  current  model,  an 
alternative  approach  using  a  moving  mesh  formulation  [2]  was  investigated.  This  would  allow 
the  mesh  to  deform  with  time  as  the  boat  moves  toward  or  away  from  a  boundary  or  other 
vehicle,  as  would  be  the  case  for  littoral  or  two-body  problems.  Further  development  is  required 
to  complete  implementation  of  the  moving  mesh  capability  in  the  6  DOF  simulation. 

The  CFX  simulation  has  been  passed  to  UNB  for  evaluation.  It  will  use  the  DRDC  generic 
submarine  shape  for  which  extensive  experimental  data  are  available.  The  maneuver  to  be 
evaluated  is  a  submarine  rising  maneuver  that  results  in  an  instability  conventional  quasi¬ 
steady  coefficient  based  hydrodynamic  models  so  far  have  not  satisfactorily  reproduced  [4,5]. 
It  is  a  good  test  case  as  the  submarine  can  be  modelled  as  an  isolated  deeply  submerged  body 
throughout  the  maneuver  (the  underwater  roll  instability  develops  before  the  boat  surfaces). 

Emergency  Rising  Maneuver 

In  emergencies,  submarines  can  blow  ballast  at  depth  to  get  to  the  surface  fast.  In  a  typical 
scenario,  ballast  is  blown  while  the  boat  is  proceeding  ahead  in  straight  and  level  flight.  Stern- 
plane  control  is  used  to  immediately  pitch  up  the  nose  of  the  boat  and  speed  is  increased  to 
minimize  large  negative  angles  of  attack  resulting  from  buoyancy.  Propulsion  and  the  buoyancy 
component  in  the  axial  direction  accelerate  the  boat  towards  the  surface  (Figure  1). 
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The  buoyancy  component  normal  to  the  hull  axis  generates  a  crossflow  (—w  in  Figure  1) 
which  results  in  large  flow  incidence  (0)  angles.  In  trials  with  small  to  medium  sized  submarines, 
this  flow  incidence  can  result  in  a  roll  instability  [4,5].  The  instability  occurs  because  the  sail 
is  pointing  into  the  crossflow.  With  small  submarines  (around  1000  t),  underwater  roll  angles 
as  large  as  25  degrees  are  seen.  With  large  diesel  boats  (2000  to  3000  t),  the  underwater 
roll  angle  is  less  than  half  that  but  still  large  enough  to  instigate  excessive  roll  when  the 
boat  surfaces  and  temporarily  loses  static  stability  (until  flood  water  drains  from  its  sail  and 
deck  casing).  The  problem  is  most  severe  in  small  boats  because  the  instability  is  caused  by 
the  destabilizing  hydrodynamic  force  on  the  sail  (proportional  to  sail  size)  overcoming  static 
stability  (proportional  to  boat  mass).  And  sail  size  tends  not  to  diminish  as  boats  get  smaller 
whereas,  of  course,  boat  mass  does. 

The  objective  of  the  evaluation  is  to  reproduce  through  simulation  the  underwater  roll 
instability  of  a  3000  t  boat  and,  if  successful,  to  investigate  mitigation  strategies.  This  report: 

1)  develops  the  solid  body  model  that  CFX  will  call  for  its  unsteady  boundary  conditions, 

2)  develops  a  fast,  quasi-steady,  coefficient  based  hydrodynamic  model  as  a  stand-in  for  the 

RANS  model  and  uses  it  to  establish  preliminary  rising  scenarios. 

The  first  step  is  to  present  the  solid-body  equations  of  motion  for  the  vehicle.  These  are 
six,  second  order,  ordinary  differential  equations  (ODEs)  adapted  from  Feldman’s  standard 
equations  [6].  The  coefficient  based  hydrodynamic  model  is  then  presented  followed  by  a  high 
incidence  propulsion  model.  A  weight  and  buoyancy  model  accounting  for  high  pitch  angles 
and  a  variable  BG  value  during  the  blow  is  developed  and,  finally,  a  control  model  is  presented 
that  is  used  for  changing  appendage  deflection  and  propulsion  states  in  the  simulations. 

The  report  finishes  by  carrying  out  several  simulations  using  the  coefficient  model.  This 
exercises  all  but  the  RANS  modelling  and  establishes  suitable  scenarios  for  the  RANS  simula¬ 
tions  to  follow.  This  is  an  invaluable  aid  in  understanding  the  complex  interplay  between  the 
many  parameters,  allowing  wise  use  to  be  made  of  the  extensive  computational  resources  the 
RANS  simulations  will  require. 

The  coefficient  based  simulations  show  a  roll  instability  is  present  and,  surprisingly,  gen¬ 
erates  significant  roll  angles  if  an  appreciable  initial  heel  angle  (~  2  degrees)  is  present.  This 
is  a  new  result  and  reflects  the  fact  that  several  key  characteristics  are,  for  the  first  time,  being 
modelled  simultaneously.  If  the  RANS  simulations  can  validate  these  results,  the  rising  stability 
problem  may  be  solved. 
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2  The  Solid  Body  Equations  of  Motion 

The  standard  submarine  equations  of  motion  [6]  are  well  known,  extensively  used  in  underwater 
vehicle  simulations,  and  readily  integrated  numerically  using  a  standard  ODE  solver.  The 
equations  describe  the  time  dependency  in  12  defining  vehicle  states: 

y  =  u,v,w,p,q,r,x0,y0,z0,<l>,6,i/)  (1) 

The  first  six  states  are  the  translational  and  rotational  vehicle  velocities  in  body  fixed  axes. 
The  next  three  are  the  coordinates  of  the  origin  of  the  body  fixed  axes  in  inertial  space.  The 
final  three  are  the  Euler  angles  orienting  the  body  fixed  axes  relative  to  inertial  space.  The 
Euler  angles  must  be  applied  in  a  consistent  order:  if  the  x,  y,  z  body  axes  are  initially  aligned 
with  the  x0l  y0,  z0  inertial  axes,  the  body  axes  are  oriented  by  yawing  about  the  z  axis  through 
an  angle  pitching  about  the  y  axis  through  an  angle  9,  and  then  rolling  about  the  x  axis 
through  an  angle  4>. 

The  six  equations  of  motion  follow.  They  are  known  as  the  ‘solid  body’  equations  in  the 
current  collaboration  as  they  describe  the  motion  of  a  rigid,  solid  body  through  a  fluid,  as 
opposed  to  the  equations  of  motion  for  the  fluid  itself  which  are  solved  by  the  RANS  solver. 
The  LHS  of  the  equations,  the  ‘rate  of  change  of  momentum’  terms  expressed  in  body  fixed 
coordinates,  are  exact  for  a  rigid  body.  The  RHS  terms  describe  the  forces  on  the  vehicle:  the 
hydrodynamic  forces  FH,  the  static  weight  W  and  buoyancy  B  forces,  and  where  applicable 
the  appendage  control  forces  Fc  and  propulsor  forces  FP. 

Axial  Force 

m  u  —  vr  +  wq  —  xG(q2  +  r 2)  +  yG{pq  —  r)  +  zG{pr  +  q ) 

=  XH(t,y)  —  (W  -  B)  sin  9  +  Xc (; t ,  y)  +  XP(t,  y)  (2a) 

Lateral  Force 

m  v  —  wp  +  ur  —  yG(r 2  +  p2)  +  zG(qr  —  p)  +  xG(qp  +  f) 

=  Ftf(t,y)  +  {W  -  B)  cos  9  sin  (f>  +  Yc(t,y)  (2b) 

Normal  Force 

m  w  -  uq  +  vp-  zG(p 2  +  q2)  +  xG(rp  -  q)  +  yG(rq  +  p) 

=  ZH(t,  y)  +  (W  -  B)  cos  9  cos  <p  +  Zc(t,y)  (2c) 

Rolling  Moment 

4 P+(Iz~Iy)(ir-(f+pq)Izx  +  (r2-q2)Iyz  +  (pr-q)Ixy+m  yG(w-uq+vp)-zG(v-wp+ur) 

=  KH(t ,  y)  +  (yGW  —  yBB )  cos  9  cos  4>  ~  {zqW  —  zbB)  cos  9  sin  (f>  +  KP(t ,  y)  (2d) 

Pitching  Moment 

IyQ+{Ix-Iz)rp-(p+qr)Ixy  +  (p2-r2)Izx  +  {qp-r)Iyz+m  zG(u-vr+wq)-xG(w-uq+vp) 

=  MH(t,  y)  —  (xGW  —  xbB)  cos  9  cos  4>  —  (zgW  —  zbB)  sin  9  +  Mc(t,  y)  (2e) 

Yawing  Moment 

hr  +  (Iy-Ix)pq-{q+rp)IyZ  +  {q2 -p2)Ixy  +  {rq-p)IZx  +  m  xG(v-wp+ur)-yG(u-vr+wq ) 

=  NH(t,  y)  +  (xGW  -  xbB)  cos  9 sin  cj)  +  (yGW  -  yBB)  sin#  +  Nc(t,  y)  (2f) 
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The  m  and  W  parameters  here  refer  to  the  total  mass/weight  enclosed  by  the  hydrodynamic 
envelope,  including  any  free  flooding  water  enclosed  by  this  envelope  and  forced  to  move  with  the 
vehicle  (eg,  ballast  tank  flood  water).  The  IVJ  terms  are  the  moments  and  products  of  inertia  of 
this  mass  and  xG,yG,  zG  locate  its  center,  all  in  body  fixed  coordinates.  The  equations  assume 
that  the  mass  m,  including  the  flood  water,  translates  and  rotates  as  a  rigid  body.  The  center 
of  buoyancy  coordinates  xB,yB,zB  locate  the  centroid  of  the  hydrodynamic  envelope.  The 
buoyancy  B  is  the  volume  within  this  envelope  multiplied  by  the  water  density. 

The  equations  of  motion  assume  the  mass  is  rigid  and  constant.  They  neglect  the  contribu¬ 
tion  of  dm/dt  and  dl/dt  terms  in  the  momentum  equation  which  occur  when  ballast  is  blown. 
This  is  justified  on  the  basis  that  the  overall  mass  change  is  small  (less  than  10%)  and  takes 
place  slowly.  The  simulations  do  model  the  change  of  mass  with  time  but  only  in  a  quasi-steady 
manner.  (The  RANS  hydrodynamic  model  is  a  true  unsteady  model.) 

Equations  (2)  are  first  order  ODE’s  in  the  body  axes  velocities  but,  implicitly,  are  second 
order  in  positional  coordinates.  Body  axis  positional  coordinates  are  not  of  interest  so  the 
following  ‘auxiliary’  first  order  ODE’s  are  integrated  simultaneously  with  (2)  to  give  the  inertial 
coordinates  and  Euler  angles: 


Xq  =  u  cos  9  cos  ip  ±  u(sin  cp  sin  9  cos  ip  —  cos  (p  sin  ip)  +  tc(sin  <p  sin  ip  +  cos  (p  sin  9  cos  VO  (3a) 

Vo  =  u  cos  9  sin  ip  +  u(cos  <p  cos  ip  +  sin  <p  sin  9  sin  VO  +  w(cos  (p  sin  9  sin  ip  —  sin  V  cos  VO  (3b) 

z0  =  —u  sin  9  +  v  cos  9  sin  (p  +  w  cos  9  cos  <p  (3c) 

4>  =  p  +  (r  cos  cp  +  q  sin  cp)  tan  9  (3d) 

9  =  q  cos  (j>  —  r  sin  <p 
■  r  cos  <p  +  q  sin  cp 

V>  =  - n - " 

cos  9 

Note  that  (3d)  and  (3f)  make  the  equations  singular  at  pitch  angles  of  ±90  degrees. 

Equations  (2)  and  (3)  define  12  nonlinear,  coupled,  first  order  ordinary  differential  equa¬ 
tions  in  the  12  states  (1).  To  integrate  the  equations  numerically,  they  must  be  in  the  form: 

y  =  f(t,  y) 


(3e) 

(3f) 


Equations  (3)  are  already  in  this  form  but  (2)  are  not.  They  are  put  in  the  correct  form  by 
inverting  the  coefficient  matrix  in  the  following  reformulation  of  (2): 
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where  ./>(t,  y)  contains  all  the  terms  from  both  the  right  and  left  hand  sides  of  (2)  that  have 
no  explicit  y  terms. 

For  simulations  in  which  the  mass,  mass  centroid,  and  moments  of  inertia  do  not  change 
with  time,  the  coefficient  matrix  in  (4)  need  only  be  inverted  once.  For  simulations  in  which 
the  mass  terms  change  with  time,  as  they  do  when  blowing  ballast,  matrix  inversion  must  occur 
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at  each  time  step.  This  is  a  straightforward  calculation  that  can  be  carried  out  quickly  using 
efficient,  compiled  algorithms  readily  available  in  many  scientific  utility  software  packages. 

In  the  RANS  simulations,  FH  in  (2)  is  modelled  in  its  entirety  by  the  RANS  solver.  This  is 
accomplished  using  a  body  fitted  mesh  about  DRDC’s  generic  submarine  shape  which  consists  of 
an  axisymmetric  hull,  a  sail,  and  four  identical  tail  fins  in  a  symmetric  ‘+’  configuration.  There 
is  no  propeller  and  so  propulsion  forces  are  accounted  for  separately.  Similarly,  the  appendages 
do  not  deflect  so  control  forces  must  be  accounted  for  separately.  Hence  the  separate  Fc  and 
FP  terms  in  (2).  Keeping  propulsion  and  appendage  deflection  out  of  the  RANS  model  will 
considerably  reduce  computation  time  and  complexity;  these  capabilities  can  be  added  in  future 
if  required.  Of  course,  the  hydrodynamic  stability  provided  by  the  appendages  is  present  in  the 
RANS  model  because  the  appendages  are  present. 

Submarine  tailfins  in  a  +  configuration  typically  do  not  deflect  differentially.  Thus,  there 
are  no  appendage  control  forces  in  the  rolling  moment  equation  of  motion. 

Using  Hydrodynamic  Coefficients  With  the  Solid  Body  Equations 

For  rapid,  preliminary  simulations,  FH  in  (2)  is  replaced  with  a  coefficient  based  model.  This 
model  accounts  for  ‘added  mass’  (see  Watt  [7])  using  acceleration  coefficients.  This  approach 
is  suggested  by  potential  flow  theory  where  the  unsteady  component  of  the  forces  exerted  on  a 
vehicle  moving  through  an  inviscid  fluid  can  be  written  exactly  as: 


Unsteady  component  of  F„  = 
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(5) 


The  coefficient  matrix  here  is  the  added  mass  matrix.  It  is  symmetric.  When  divided  by  p  it 
is  a  function  only  of  vehicle  geometry.  Potential  flow  predictions  of  these  coefficients  agree  well 
with  measurements  made  in  simple  acceleration  (oscillation)  experiments  in  water.  On  the  other 
hand,  potential  flow  predictions  of  the  steady  state  forces  do  not  agree  well  with  experiment 
because  they  cannot  predict  vorticity  nor,  therefore,  lift  on  a  moving  body.  Vortical  flows  also 
result  in  unsteady  forces  that  are  not  modelled  by  (5)  and  which  may  be  a  contributing  factor 
to  the  failure  of  coefficient  based  models  to  predict  the  rising  instability  [5]. 

The  diagonal  added  mass  matrix  coefficients  have  the  largest  magnitudes.  The  following 
off-diagonal  coefficients  are  identically  zero  when  the  vehicle  has  a  vertical  plane  of  symmetry 
as  the  current  geometry  does: 


XXX 


Y  Y  Y 

±  U  5  W  5  J  q  1 


ZrinZ^Zf.,  Ku,  Kw,  K;. .  M^,  M„, Mj.,  N^N^N^O 


(6) 


The  reason  for  introducing  the  coefficient  model  in  this  section  is  that  its  added  mass  model 
necessitates  changes  to  (4).  Explicit  acceleration  terms  in  FH  in  (2)  must  be  combined  with 


DRDC  Atlantic  TR  2007-008 


5 


similar  terms  on  the  LHS  so  that  (2)  becomes: 
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(7) 


As  for  (4),  the  coefficient  matrix  in  (7)  must  be  inverted  at  each  timestep  of  the  numerical 
integration  if  the  mass  terms  change.  The  added  mass  coefficients  do  not  change  as  long  as  the 
external  vehicle  shape  does  not  change. 


3  The  Quasi- Steady,  Coefficient  Based  Hydrodynamic  Model 

This  model  is  used  for  preliminary  calculations,  to  explore  the  effects  of  the  modelling  pa¬ 
rameters  on  the  simulation  and  select  values  which  reproduce  maneuvers  with  the  desired 
characteristics.  It  is  based  on  Feldman’s  standard  equations  [6].  Differences  are: 

•  except  for  added  mass,  Feldman’s  unsteady  model  is  not  used, 

•  translation  hydrodynamics  are  modelled  using  high  incidence  experimental  data  acquired 
at  up  to  30  degrees  incidence  at  R  =  23  million  with  the  same  generic  submarine  shape 
used  in  the  RANS  model;  this  data  is  presented  in  Figure  8  in  Reference  [5]  and  referred 
to  below  as  the  Fuvw  function, 

•  a  high  incidence  propulsion  model  is  used  based,  again,  on  experiments  with  our  generic 
submarine  shape, 

•  simple  tailplane  control  models  are  used  which  allow  the  desired  maneuver  to  be  achieved; 
the  effect  of  propulsive  state  on  tailplane  control  is  not  modelled. 

The  subsections  below  present  the  gF  functions  from  (7),  isolating  the  quasi-steady  terms 
from  the  LHS  of  (7)  on  the  first  line,  and  the  FG  functions  from  the  RHS  of  (2)  so  they  are 
easily  identified  for  use  in  the  RANS  simulation.  The  propulsion  and  weight  and  buoyancy 
models  are  discussed  in  detail  in  subsequent  sections. 

The  last  lines  in  the  gF  functions  below  contain  special  functions  and  the  weight  and 
buoyancy  terms.  The  middle  lines  contain  constant  coefficients  multiplying  state  velocities. 
The  coefficients  are  either  suggested  by  Feldman,  by  potential  flow  analysis  [7],  and/or  result 
from  force  estimates  made  by  Mackay’s  DSSP20  program  [8]  and  converted  to  first  and  second 
order  derivatives  using  the  DERIVS  program  [9] . 

The  Fuvw  functions  are  plotted  in  Figure  2.  These  functions,  all  of  the  coefficients  listed 
below,  the  mass  and  added  mass  coefficients  from  the  LHS  of  (7),  and  many  other  relevant 
physical  constants  for  the  simulations  are  listed  in  Appendix  A.  All  of  the  F%3  coefficients 
below  are  constants  unless  a  dependency  is  explicitly  indicated. 
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Axial  Force 

Feldman’s  Xvvv2  +  Xwww2  —  Drag  translation  terms  are  included  in  the  X'uvw  ( 0 ,  <J>)  function. 
To  avoid  having  to  model  serious  discontinuities,  this  function  ignores  the  effect  of  sail  stall 
which  occurs  when  0  >  23  degrees  for  <h  =  60  to  120  degrees,  a  flow  regime  not  expected  to 
be  important  for  the  current  simulations.  The  only  other  force  suffering  the  same  restriction  is 
rolling  moment. 

9x(t,  y)  =  m[vr  -  wq  +  xG(q2  +  r2)  -  yGpq  -  zGpr\ 

+  Xuquq  +  Xvrvr  +  Xwpwp  +  Xwqwq  +  Xppp2  +  Xrprp  +  Xqqq2  +  Xrrr 2  +  XqMq\q\ 

+  ^pU2i2X'uvw(e^)  +  Xc{t,y)  +  Xp(t,y)  -  (W-B)  sin  6  (8a) 

where: 

Xc(t,  y)  =  (X'5sSads2  +  X'Sr&r52)  u 2  (8b) 

These  latter  terms  model  the  drag  generated  by  tailplane  deflections. 

Lateral  Force 

A  linear  model  is  used  for  tailplane  control  derivatives  which  ignores  tailplane  stall.  The  lateral 
force  generated  by  the  rudder  is  not  large  but,  since  the  rudder  is  at  the  end  of  the  boat,  the 
long  moment  arm  provides  good  yaw  control. 

9v(t,y)  =  m[wp  —  ur  +  yG(r 2  +  p2)  -  zGqr  -  xGqp\ 

Yupup  +  Yurur  +  Ywpwp  +  Ywrwr  +  Ypqpq  +  Yqrqr  +  YpMp\p\  +  Yr\r\r\r\ 

+  ^pU2£2Y:vw(Q,<S>)  +  Yc(t,y)  +  (W  -  £?)  cos  0  sin  </>  (9a) 

where: 

Yc(ty)  =  -2p£2Ylu2Sr  (9b) 

Normal  Force 

Sternplane  control  is  to  normal  force  what  rudder  control  is  to  lateral  force.  It  is  a  small  force 
that  takes  advantage  of  a  long  moment  arm  to  provide  good  pitch  control. 

y)  =m[uq-vp  +  zG(p 2  +  q2)  -  xGrp  -  yGrq] 

^vp^P  “I-  Zwpwp  +  zwqwq  +  ZppP  +  Zrprp  -\-  Zqqq  -1-  Zrrv  -\-  Zq^q^q\q\ 

+  1 PU2eZ'uvw{Q ,  T)  +  Zc(t,  y)  +  (W  -  B)  cos  9  cos  <f>  (10a) 

where: 

Zc{t,y)  =  ±pPZ'su26a  (10b) 
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Figure  2  The  six  Fuvw  (0,  <I>)  forces  used  in  the  coefficient  model  (angles  in  degrees). 
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Rolling  Moment 

The  K's  term  is  included  for  generality  but  it  is  zero  for  our  generic  shape  which  has  a 
symmetric  rudder  centered  on  a  symmetric  hull  and  a  body  axes  origin  located  on  the  hull 
centerline.  The  propulsion  system  uses  a  propeller  which  uses  torque  to  generate  thrust.  This 
generates  a  moment  KP  capable  of  rolling  the  boat  a  fraction  of  a  degree  or  so.  It  is  useful 
to  model  in  a  rising  stability  investigation  because  it  introduces  a  small  roll  angle,  the  initial 
displacement  a  roll  instability  needs  to  take  hold.  As  with  axial  force,  the  K!avw(@,  d>)  function 
does  not  model  the  discontinuity  resulting  from  sail  stall  which  occurs  when  0  >  23  degrees 
for  <1>  =  60  to  120  degrees. 

Heel  and  stability  are  controlled  by  the  lateral  and  vertical  locations  of  the  centers  of 
gravity  and  buoyancy.  In  general,  yB  «  0  and  is  fixed  by  the  shape  of  the  hydrodynamic 
envelope.  On  the  other  hand,  yG  can  be  adjusted  by  moving  mass  laterally  which  is  typically 
done  to  trim  the  boat  in  roll.  The  important  ( zbB  —  zGW)  cos#sin<)>  static  stability  term  is 
the  only  mechanism  the  boat  has  for  remaining  upright  in  the  presence  of  propeller  torque  and 
destabilizing  hydrodynamic  forces.  Static  stability  is  proportional  to  BG  =  zG  —  zB,  the  height 
of  the  CB  above  the  CG. 

gx(t,  y)  =  (4  -  4 )qr  +  hxPQ  -  4 z(r2  “  #2)  ~  4 yW  +  m[yG(uq  -  vp )  -  zG(wp  -  ur )] 

+  Kupup  +  Kurur  +  Kvqvq  +  Kwpwp  +  Kwrwr  +  Kpqpq  +  Kqrqr  +  KpMp\p\  +  Kr\rf\r\ 

+  -pH2 £3 K'uvw(@ ,  d>)  +  Kc(t,y)  +  RP(t,  y) 

+  ( yGW  —  yBB )  cos  0  cos  4>  —  ( zGW  —  zbB)  cos  9  sin  0  (Ha) 

where: 

Kc(t,y)  =  ^p£3K'Sru25r  (lib) 


Pitching  Moment 

The  important  M$g  coefficient  defines  linear  pitch  control  and  allows  for  second  order 

effects  in  lift  which  can  occur  in  low  aspect  ratio  appendages;  this  latter  term  does  not  model 
stall.  The  M$3gs ,  Mgr$r  coefficients  account  for  any  moment  generated  by  the  drag  terms  in 
(8b)  not  acting  through  the  body  axes  origin,  so  they  are  zero  for  the  generic  model  with  its 
symmetric  tailplanes.  Again,  the  (zbB  —  zGW)  sin#  static  stability  term  tries  to  keep  the  boat 
level  in  pitch.  Pitch  trim  is  controlled  dynamically  using  5S  or  statically  by  adjusting  the  axial 
location  of  the  CG. 


5m(*, y)  =  (4  -  4 )rp  +  Ixyqr  -  Izx{p2  -  r 2)  -  Iyzqp  +  m[zG(vr  -  wq )  -  xG(uq  -  vp)] 

+  Muquq  +  Mvpvp  +  Mvrvr  +  Mwpwp  +  Mwqwq  +  Mppp 2  +  Mprpr  +  Mqqq2  +  Mrrr2  +  Mq^q\q\ 

+  -pU2 13 M'uvw(® ,  d>)  +  Mc(t,  y)  —  {xGW  —  xbB )  cos  9  cos  0  —  (zGW  —  zbB)  sin#  (12a) 


where: 


Mc(t,  y)  —  -pi3  ^MgsSa  +  M'slSsl8s\6s\  +  M'5s5s82  +  M'5r5rSr2^j 


(12b) 
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Yawing  Moment 

The  weight  and  buoyancy  terms  in  yawing  moment  are  zero  for  a  neutrally  buoyant,  well 
trimmed  boat.  They  do  not  interact  with  yaw  angle  or  yaw  control. 

Figure  2  shows  that  pitching  and  yawing  moments  both  have  considerable  variation  at 
high  incidence.  Since  these  functions  result  from  least  squares  fits  to  data  taken  in  <J>  increments 
of  30  degrees,  this  variation  may  not  be  adequately  resolved. 

9N(t,  y)  =  (4  -  Iy)pq  +  Iyzrp  -  Ixy(q 2  -  p2)  -  Izxrq  +  m[xG(wp  -  ur)  -  yG(vr  -  wq)\ 

+  Nupup  +  Nurur  +  Nvqvq  +  Nwpwp  +  Npqpq  +  Nqrqr  +  NpMp\p\  +  Nr]r\r\r\ 

+  \pU2tiN'uvw(®,§)  +  Nc(t,  y)  +  ( xGW  -  xbB)  cos  6  sin  <j>  +  (; yGW  -  yBB)s\n6  (13a) 
where: 

Nc(t,y)  =  ^pl3N'Sru28r  (13b) 


4  Propeller  Thrust  and  Torque 

A  high  incidence  (up  to  30  degrees)  propulsion  model  is  developed  following  Watt  [10]  by 
adapting  the  classic  propulsion  model  [ll]  so  that  wake  fraction  varies  with  incidence. 

A  conventional  dimensionless  representation  for  open  water  propeller  thrust  T  and  torque 
Q  uses  thrust  and  torque  coefficients  which  depend  solely  on  the  advance  ratio  J : 


kt(J) 


T 

pn2DA  ’ 


Kq(J) 


Q 

pn2D5  ’ 


(14) 


Va  is  the  speed  of  advance  of  the  propeller  through  the  water,  n  is  propeller  revolutions  per 
second,  and  D  is  propeller  diameter.  A  stern  propeller  on  a  submarine  operates  in  a  wake  which 
reduces  the  average  inflow  to  the  propeller  relative  to  the  speed  of  the  boat.  This  is  accounted 
for  classically  with  a  one-dimensional  correction.  The  speed  of  advance  of  the  propeller  through 
the  water  is  approximated  by: 

Va  =  (1  —  iuT)u  (15) 

where  wB  is  the  Taylor  wake  fraction  and  u  is,  as  usual,  the  forward  speed  of  the  boat.  The 
classical  model  does  not  account  for  crossflow  and  does  not  distinguish  between  u  and  U  so  (15) 
is  a  minor  extension  to  the  classical  approach.  As  shown  below,  a  further  necessary  extension  is 
to  allow  the  wake  fraction  to  vary  with  incidence,  whereas  it  is  constant  in  the  classical  model. 

In  addition,  when  a  propeller  generating  thrust  T  operates  behind  a  hull,  it  also  induces 
negative  pressure  on  the  hull  afterbody  upstream  of  it.  This  increases  the  drag  on  the  hull, 
negating  some  propeller  thrust.  This  is  accounted  for  classically  using  the  thrust  deduction 
fraction  t,  another  constant: 

XP  =  (1  -  t)T  (16) 

The  complimentary  fraction  1  —  i  is  called  the  thrust  deduction  factor. 

To  use  (14),  open  water  experiments  with  a  propeller  operating  at  various  advance  ratios 
are  required  to  determine  KT  and  Kq.  To  use  (15),  a  wake  survey  needs  to  be  conducted 
behind  the  hull  on  which  the  propeller  is  to  be  used.  And  (16)  requires  model  tests  with  the 
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Figure  3  The  DRDC  Static  Test  Rig  (STR)  generic  hull,  sail,  and  tail  with  propulsion  system 
installed  in  the  Institute  for  Aerospace  Research  9  m  wind  tunnel.  The  conventional  three  bladed 
propeller  is  mounted  on  a  six  component  internal  balance.  The  model  is  yawed  and  pitched 
using  a  floor  turntable  with  the  sail  located  vertically  or  horizontally  on  the  axisymmetric  hull. 


hull  and  propeller  combination.  Herein,  a  single  model  setup  using  the  generic  hull  with  a 
propeller  operating  behind  it  is  used  instead.  This  is  possible  only  because  of  the  information 
obtained  from  testing  at  incidence  and  because  of  the  assumptions  made  in  the  modelling  of 
incidence,  as  described  below. 

The  experiments  again  take  place  in  the  wind  tunnel  (Figure  3)  at  incidence  angles  from 
—30  to  30  degrees  at  Reynolds  numbers  over  20  million  (based  on  hull  length)  [12,13,14],  The 
propeller  diameter  is  half  the  hull  diameter  and  operates  in  dynamic  similarity  to  an  imagined 
full  scale  prototype.  Propeller  thrust  and  torque  are  measured  as  well  as  the  overall  forces  on 
the  vehicle. 

Because  and  Wj-  are  unknown,  the  effect  of  propulsion  on  KT  and  Kq  is  initially 
found  as  a  function  of: 

=  ^  =  <17> 

where  J b  is  the  behind-the-boat  advance  ratio.  The  coefficients  are  obtained  by  measuring 
propeller  thrust  and  torque  at  zero  incidence,  fitting  them  with  polynomial  interpolants,  and 
correcting  for  minor  compressibility  effects  due  to  high  propeller  RPM  [13]: 

KT(Jb )  =  0.410758  -  0.115654,4  -  0.107836Jb2  +  0.0713369Jb3  -  0.00620451  Jb4 

-  0.0127538Jb5  +  0.004878934®  -  0.000678484 J7b  +  0.00003334634®  (18a) 

KQ(.Jb)  =  0.0690631  -  0.02496584  -  0.00623472 +  0.0017180743  +  0.0057916944 

-  0.005596304s  +  0.0017895046  -  0.00024688647  +  0.00001260294®  (18b) 

These  are  plotted  in  Figure  4. 
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Figure  4  STR  zero  incidence  thrust  and  torque  coefficients  as  a  function  of  behind-the-boat 
advance  ratio.  Corrected  for  compressibility  effects.  R  >  20  million. 

Effect  of  Incidence 

Consider  Figure  5  where  STR  thrust  and  torque  data  from  two  pitch  runs  at  constant  propeller 
RPM  are  plotted.  Two  data  points  are  seen  at  each  angle  for  each  run  for  each  of  thrust  and 
torque  because  each  run  is  a  sweep  from  a  =  —30  to  30  degrees.  The  ±a  pairs  are  fairly  close 
indicating  symmetry  about  a  =  0.  Yaw  sweep  data  lose  some  of  this  symmetry  due  to,  it  is 
thought,  the  trailing  vortex  field  from  the  sail;  for  rising  stability,  vertical  plane  incidence  is  of 
primary  interest  and  so  only  the  pitch  data  are  used. 

The  wind  tunnel  speed  is  constant  for  these  pitch  sweeps.  Therefore,  as  incidence  increases, 
u  and  Jb  reduce  as  the  cosine  of  the  incidence  angle,  as  shown  in  Figure  5.  The  solid  line  fit 
through  the  Jb  data  is  the  curve: 

Jbfit(0)  =  1.02266  cos  0  (19) 

This  cos  0  variation  in  Jb  will  have  negligible  effect  on  thrust  and  torque  at  low  incidence. 
However,  in  Figure  5,  thrust  and  torque  initially  decrease  relatively  fast  with  incidence  as  it 
increases  from  zero.  This  is  thought  to  result  from  the  crossflow  sweeping  the  wake  aside  thereby 
reducing  the  wake  fraction  and  increasing  propeller  inflow.  This  reduces  propeller  loads.  The 
load  reduction  is  curtailed  at  15  to  20  degrees  incidence  because  of  the  increasingly  strong 
roll-off  in 

Thus,  Kt  and  Kq  vary  with  0  independently  of  Jb,  which  means  the  description  in  (18) 
is  incomplete.  From  basics,  it  is  known  that  thrust  and  torque  can  be  described  solely  in  terms 
of  the  conventional  advance  ratio  J  =  V^/nD  and,  therefore,  that  the  functions: 

Kt,q  ( T - ")  (20) 

\1  —  wToJ 

are  always  correct.  Here,  wTo  is  wT  at  zero  incidence  where  KTiKq  (18)  were  acquired;  it 
is  the  conventional  constant  Taylor  wake  fraction.  To  use  (20),  and  therefore  w?  must  be 
known  as  a  function  of  incidence.  Then  (17)  allows  the  thrust  and  torque  coefficients  to  be 
generalized: 

Kr,Q(Jb)  — 5 "  Kt,q  ( 7 - —  Jb)  (21) 

\  l  —  wTo  j 

In  other  words,  (17)  remains  valid  in  the  generalized  model. 
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Figure  5  Propeller  thrust,  torque,  and  behind-the-boat  advance  ratio  variation  with  incidence 
for  <P  =  0, 180  degrees,  constant  wind  speed,  and  constant  propeller  RPM  nominally  set  for  self 
propulsion.  Not  corrected  for  compressibility  effects,  o  =  Run  208  and  +  =  Run  210  from  [12]. 


as: 


The  shapes  of  the  curves  in  Figure  5  suggest  the  drop  in  wake  fraction  might  be  modelled 


VJ'j'  —  Wj'qC 


-(key 


(22) 


where  Wp0,  k,  7  are  unknown  constants.  Using  (19)  in  a  least  squares  fit  of  (21)  and  (22)  to  the 
Kt,Kq  incidence  sweep  data  (multiplying  KT  by  1.063  and  Kq  by  1.048  because  the  data 
are  uncorrected  for  compressibility  effects)  gives: 


wT  0  =  0.31,  k  =  3.4,  and  7  =  1.18 


(23) 


The  lines  through  the  Kt,Kq  data  in  Figure  5  show  the  fit. 

To  summarize,  the  effect  of  propulsion  on  axial  force  and  rolling  moment  is  modelled  using: 


XP  =  pn2D\l-t)KT(Jm) 
KP  =  -pn2D5KQ(Jm) 


where  Jm  is  the  modified  advance  ratio: 


u(  1  —  wToe  (k0P'j 
nD{l  —  wTo) 


and  the  I\t,Kq  functions  in  (24a)  are  given  by  (18). 


(24a) 


(24b) 
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a 

b)  Thrust  deduction  factor. 


Figure  6  The  thrust  deduction  data  in  (b)  are  the  differences  between  the  fitted  curves  from 
(a)  normalized  by  the  thrust  data  from  Figure  5.  o  =  Run  208,  +  =  Run  210,  •  =  Run  214, 
X  =  Run  216  from  [12]. 


Thrust  Deduction 

The  thrust  deduction  factor  1— t  =  XP/T  is  estimated  from  the  STR  experiments  by  comparing 
overall  axial  force  and  localized  propeller  thrust  measurements.  Runs  208  and  210  from  [12] 
give  overall  axial  force  X  and  propeller  thrust  T  measurements  for  pitch  sweeps  from  a  =  —30 
to  30  degrees  with  the  propeller  providing  self-propulsion  at  zero  incidence.  These  are  the  red 
data  points  in  Figures  5  and  6a.  To  get  the  XP  component  from  the  X  data,  the  axial  force 
from  similar  data  acquired  with  the  propeller  unpowered  and  windmilling  is  subtracted  from 
the  X  data  (propeller  thrust  for  the  windmilling  condition  is  effectively  zero).  The  windmilling 
X  data  (runs  214  and  216  from  [12])  are  shown  in  blue  in  Figure  6a. 

Unlike  the  yaw  data  presented  in  [14],  the  Figure  6a  pitch  data  are  uncorrected  for  support 
strut  tare  and  interference  effects.  This  correction  would  put  the  self-propulsion  zero  incidence 
axial  force  much  closer  to  zero,  as  would  be  expected.  However,  it  is  only  the  difference  between 
the  self-propulsion  and  windmilling  X  data  that  is  used  and  this  difference  is  unaffected  by  the 
correction  which  is  the  same  for  each  data  set. 


The  STR  axial  force  data  clearly  suffer  from  some  systematic  error  as  the  model  is  pitched 
through  zero  degrees.  To  correct  for  this,  and  because  the  data  in  the  two  sets  of  runs  are  not 
all  at  the  same  angles,  the  data  are  least  square  fitted  with  4th  order  polynomials  in  a  (even 
powers  only).  These  fits  are  shown  in  Figure  6a.  The  thrust  deduction  factor  data  is  then 
calculated  from: 


Xfit  -  Xfit 

208,210  nt214,216 


pn2DAKT 


208,210 


(25) 


and  is  shown  in  Figure  6b. 

Although  1  —  t  appears  to  vary  with  incidence,  the  variation  is  similar  in  magnitude  to 
the  fitting  error  in  Figure  6a  and  is  discounted.  Thus,  thrust  deduction  is  just  a  constant,  the 
average  of  the  Figure  6b  data: 

1  -t  =  0.8662  (26) 
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Commanded  Speed  and  Trim 

It  is  usually  more  convenient  to  command  forward  speed  than  RPM.  Commanded  forward  speed 
uc  is  the  self-propulsion  speed  of  the  boat,  the  speed  achieved  when  in  ‘equilibrium’  (straight, 
level,  neutrally  buoyant,  zero  incidence,  fully  trimmed  flight).  The  equilibrium  hydrodynamic 
state  is  the  same  at  any  speed  and  associated  with  it  is  a  simple  linear  relationship  between 
speed  and  RPM,  as  given  by  the  advance  ratio  (17).  To  use  this  relationship,  the  self-propulsion 
behind-the-boat  advance  ratio  Jbs  must  be  known. 

Equilibrium  implies  roll  and  pitch  angles,  flow  incidence,  all  accelerations,  and  all  velocities 
are  zero  except  for  the  self-propulsion  (commanded)  speed: 

u  =  U  =  uc  (27) 


When  equilibrium  is  achieved,  only  the  four  X,  Z,  K,  M  equations  of  motion  have  nonzero  terms 
and  these  must  combine  in  such  a  way  as  to  zero  each  equation: 


X  : 

0  =  ^pu2J2  (. X'uvw(0 ,  T)  +  X5s5Js2)  +  pn2D\  1  -  t)KT(Jbs) 

(28a) 

Z  : 

0  =  ^pu2J2  (Z'uvw(  0,  4»)  +  Z'S'6a)  +W-B 

(28b) 

K  : 

0  =  pn2D5KQ(Jbs )  +  ( yGW  -  yBB) 

(28c) 

M  : 

0  =  X-pu\e  {M'uvw (0,  4>)  +  M'Ss5s)  -  (xGW  -  xbB) 

(28d) 

(28d)  reflects  the  fact  that  M's  ^  |  =  MSsSs  =  0,  as  given  in  Appendix  A.  The  expressions  for 
Fuvw  in  Appendix  A  show  that  all  the  horizontal  plane  Y,  K,  N  hydrodynamic  translational 
forces  are  zero  when  0  =  0.  Thus,  simply  zeroing  the  rudder  deflection  zeros  the  Y  and  N 
forces.  On  the  other  hand,  the  vertical  plane  X,  Z,  M  translational  forces  reduce  to  nonzero 
constants  at  0  =  0  because  of  the  hydrodynamic  asymmetry  associated  with  the  sail. 

Equations  (28)  show  that  trimming  is  necessary  to  achieve  equilibrium.  The  asymmetric 
hydrodynamic  forces  are  normally  countered  using  hydrodynamic  forces  generated  by  the  stern 
and  bowplanes;  this  keeps  the  trim  effective  as  speed  is  changed.  However,  since  bowplanes  are 
not  present  in  the  current  model,  trim  is  achieved  instead  using  a  combination  of  sternplane 
deflection  and  weight  compensation.  A  small  sternplane  deflection  is  used  to  trim  the  boat  to 
zero  pitch  angle  and  the  weight  is  adjusted  to  negate  the  normal  force  from  the  sternplanes. 
The  weight  trim  Wt  is  added  at  xGo  =  xB  =  0  so  that  the  last  term  in  (28d)  remains  zero. 
The  sternplane  trim  deflection  is  therefore  independent  of  speed: 


$st  = 


-M'uvw{  0,$>)  180 


Mi 


7 r 


=  0.813  degrees. 


(29a) 


while,  from  (28b),  the  weight  trim  is  speed  dependent: 


Wt  =  W  -  B  =  pu2cl 2  (Z'uvw (0,  $)  +  Z'&s5st )  =  -543 u2c  Newtons 


(296) 


At  5  knots,  Wt  is  just  0.01%  of  B  so  its  speed  dependence  minimally  compromises  equilibrium. 

The  sternplane  trim  deflection  is  so  small  it  has  only  a  small  effect  on  the  solution  of  (28a) 
for  Jbs : 

Jbs  =  0.9783  (30) 
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With  the  unique  constant  Jbs  known,  the  desired  relationship  between  n  and  uc  is: 


n  = 


uc 

JbsD 


0.2555 uc 


(31) 


and  (24b)  can  be  replaced  with: 


^ m 


uJhs(l-wToe  (fcQ)7) 
W(1  -  Wto) 


(32) 


In  other  words,  the  advance  ratio  is  just  proportional  to  the  actual  forward  speed  divided  by 
the  forward  speed  that  the  current  propeller  RPM  achieves  in  equilibrium  conditions. 

Finally,  consider  (28c)  which  states  that  propeller  torque  must  be  trimmed  with  a  lateral 
shift  in  the  CG  to  achieve  equilibrium.  This  will  not  be  done  in  the  current  simulations  so 
that  a  small  roll  (heel)  angle  will  always  be  present,  as  is  typical  in  real  life.  This  should  allow 
any  roll  instability  that  is  present  to  develop  and  grow.  If  yG  =  yB.  then  propeller  torque  will 
generate  a  small  roll  angle  and  (28c)  is  more  properly  modelled  as: 


K  : 


0  =  pn2D5KQ(Jbs)  -  (. zGW  -  zbB)  sin  <j> 


(33) 


which  gives: 


<f>0  =  —0.1044  degrees. 


at  a  steady  forward  speed  of  5  knots.  This  increases  with  n2  and  therefore  with  u 


2 

C  * 


5  Buoyancy,  Weight,  and  Blowing 

Buoyancy  is  fixed  by  the  shape  of  the  external  hydrodynamic  envelope.  With  the  main  ballast 
tanks  (MBTs)  flooded,  the  weight  is  nominally  equal  to  the  buoyancy.  However,  as  noted 
above,  the  weight  can  be  adjusted  (using  small  onboard  compensation  tanks)  for  each  steady 
state  condition  at  which  the  operators  wish  to  achieve  equilibrium.  The  simulation  initial 
conditions  account  for  this  trim  so  the  simulations  begin  in  equilibrium.  Blowing  must  also  be 
modelled,  but  it  is  always  activated  after  a  simulation  has  been  started  in  equilibrium. 

Weight  is  modelled  as: 


W  =  W0-[iB 

where 

W0  =  B  +  Wt 

Similarly: 

II 

O 

and 

B 

m0  = - lmt 

9 

(34a) 

(34b) 


The  ‘o’  subscript  refers  to  conditions  immediately  prior  to  blowing.  Wa  is  the  equilibrium 
weight  which  equals  the  buoyancy  plus  the  trim  weight  Wt  necessary  to  achieve  equilibrium  at 
some  initial  state.  The  blown  mass  fraction  y  is  zero  at  submerged  equilibrium  and  about  0.1 
with  the  MBTs  empty.  As  was  shown  above,  the  trim  weight  Wt  is  approximately  0.0001.B  so 
Wa  is  effectively  equal  to  B.  Since  rotation  is  zero  during  equilibrium,  the  values  of  the  moments 
of  inertia  are  not  critical  and  any  effect  Wt  has  on  the  moments  is  ignored.  It  is  assumed  that 
both  W0  and  Wt  have  their  centroids  at  xGo,yGo,  zGo  and  that  initial  equilibrium  in  straight 
and  level  flight  is  achieved  with  xGo  =  xB  and  yGo  =  yB. 
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Mass  Change  While  Blowing 

Submarines  have  several  main  ballast  tanks  distributed  along  the  length  of  the  hull  and  can  blow 
them  independently  or  all  together.  The  forward  tanks  are  always  largest  so  the  nose  rises  faster 
than  the  stern.  This,  coupled  with  propulsion,  gets  the  boat  to  the  surface  fast  in  an  emergency. 
Emergency  operating  procedures  may  call  for  the  forward  tanks  to  be  blown  first  and  the  aft 
tanks  only  after  the  nose  has  begun  to  rise.  However,  all  tanks  are  blown  simultaneously  in 
the  current  model.  Each  tank  centroid  is  assumed  to  be  located  on  the  hull  centerline.  Since 
the  forward  MBT  volumes  are  largest,  the  blown  mass  centroid,  (x^,  y^,  ),  will  be  slightly 

forward  of  the  CB.  There  are  N  main  ballast  tanks  with  a  total  volume  V B  =  Y^iLi  V with 
each  tank  having  an  axial  centroid  at  x Bi.  The  blown  mass  fraction  is  then: 


P  = 


WD  -  W 
B 


N 

pg  v«i 

i= 1 


B 


N 

i= 1 


(35) 


where  Va  is  the  total  volume  of  air  in  (water  expelled  from)  the  MBTs  and  Vai  is  the  volume  of 
air  in  tank  i.  When  blowing  begins,  Vai  =  0  and  it  is  maximum  when  all  the  water  is  expelled, 
when  Vai  =  VT%. 

The  axial  centroid  of  the  blown  mass  fraction  is: 


N  N 

XT  iPi  XTiVai 


= 


i=l 


i= 1 


P 


vn 


(36) 


However,  the  vertical  centroid  of  the  blown  mass  fraction  is  dependent  on  the  local  vertical 
centroids,  the  z^.  say,  which  vary  approximately  linearly  with  local  blown  mass  fraction 
even  though  the  MBTs  have  irregular  shapes.  The  z^.  have  an  initial  value  near  the  top  of 
the  hull  when  //,;  =  0  and  end  with  a  value  on  the  centerline  when  the  blow  is  complete.  The 
initial  zfl/.  value  is  taken  to  be  90%  of  the  maximum  hull  radius  d/2: 


zn  = 


N 

i= 1 _ 

p 


N 


-0.45dJ^  (  1 


i= 1 


Va- 

Vt 


Vn 


Vn 


(37) 


This  accounts  for  different  levels  of  water  in  each  MBT,  which  occurs  when  the  boat  is  pitched 
up  putting  the  forward  tanks  at  a  different  depth,  and  therefore  pressure  and  blown  air  density, 
than  the  aft  tanks. 

Expressions  for  xG,  zG  in  terms  of  z M  are  obtained  by  taking  moments  about  the  body 
axes  origin.  Remembering  that  B  and  xB  do  not  change  when  blowing  ballast  and  xGo  =  xB, 
zGo  =  zB  +  BGa  initially,  there  results: 


xG  = 


jiB 

x B  ~  WnXfl 


1  - 


fj,B 

wn 


ZG  ~ 


zB  +  BG0  -  ^~zfj 


1  - 


l^iB 
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(38) 
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(39) 


The  BG  =  zG  —  zB  value  is  then: 


BG 


BG0  +  zB  -  z^) 


pB 


pB 

BG0  ~l  ~fiB^ZGo  ~  Z ^ 
1_  Wo 


If  the  vehicle  CG  is  on  the  hull  centerline,  and  because  the  vertical  centroids  of  the  MBTs  are 
on  the  hull  centerline,  then  BG  has  the  same  value  before  (when  p  =  0)  and  after  the  blow 
(when  z p  =  0).  During  the  blow,  z p  is  negative  so  BG  temporarily  increases. 

From  (2d),  static  stability  in  roll  is  determined  by  the  term  zGW  —  zbB  which,  using  (34a) 
and  (38b),  becomes: 

ZqW  -  zbB  =  B  ( BG0  -  pZp)  =  B  BG*  (40a) 

where: 

BG*  =  BG0  -  pZfl  (406) 

When  W  ^  B,  zGW  —  zbB  is  different  from  B  BG  or  WBG  and  it  is  convenient  to  use  the 
modified  BG  value  BG  .  BG  is  the  same  before  and  after  the  blow,  regardless  of  zGo,  and 
differs  from  BG  only  by  0(/xzGo ,  p2 z M).  The  temporary  increase  in  BG  during  the  blow  tends 
to  delay  the  onset  of  the  roll  instability  until  the  MBTs  are  almost  empty.  So  modelling  the 
vertical  variation  of  z;j  with  time  is  important. 

The  moments  of  inertia  all  vary  somewhat  with  blown  mass.  These  variations  are  small 
and  often  neglected  during  blowing  simulations,  primarily  so  the  mass  matrix  multiplying  the 
6  accelerations  in  the  equations  of  motion  does  not  have  to  be  continually  inverted  during 
numerical  integration  of  the  equations.  However,  the  variation  is  easy  to  account  for  and 
inverting  a  6  x  6  matrix  is  not  computationally  expensive  these  days.  Here,  then,  are  the 
moments  and  products  of  inertia  accounting  for  blown  mass  by  assuming  it  occurs  at  a  point: 


4 


4 


J  ( V 2  +  z2)dm  =  Ixo 
J  (z2  +  x2)dm  =  I y  0 
[  (x2  +  y2)dm  =  Izo 


/xm0z2 

ixm0(z2  +  xl) 

ixm0xl 


IXy  =  I  xy  dm  =  I 


xy0 


IVz  =  yz  dm  =  / 


yzo 


(41) 


4x  =  zx  dm  =  Izxo  +  nm0z„x 


M  M 


Blowing  Model 

The  above  mass  model  requires  knowledge  of  Vai,  the  volume  of  air  in  each  MBT.  This  volume 
is  derived  here  as  a  function  of  time,  tank  location,  and  the  depth  and  pitch  angle  of  the  boat. 


The  MBTs  are  blown  from  a  reservoir  of  very  high  pressure  air  (several  hundred  atmo¬ 
spheres)  which  discharges  through  nozzles  in  the  MBTs.  One  dimensional,  isentropic,  compress¬ 
ible  flow  theory  [15]  predicts  the  maximum  nozzle  velocity  to  be  Mach  1.  This  is  maintained 
if  the  discharge  to  reservoir  pressure  ratio  is  less  than  0.53  which  will  be  the  case  for  most  of 
the  duration  of  any  blow.  In  this  case,  the  reservoir  air  density  pr  during  the  blow,  assuming 
isentropic  flow,  is: 


dpr(t) 

dt 


(42) 
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where  k  =  1.4  for  air  and  is  a  negative  constant  that  depends  only  on  the  properties  of  the 
bottled  air  and  nozzle  diameter: 

dpr  I 


<?!  = 


dt 


fe+i 


(43) 


Integrating  (42)  gives  the  isentropic  blowing  model: 


pr{t) 


Pro 


k-  1 
2 


C^t 


where  C2  is  the  blowing  constant: 


fe-i 


dpr 


dt 

Pro 


o 


dmr 


dt 


mr 


(44) 


(45) 


mr  is  the  mass  of  air  in  the  reservoir  and  ta  is  the  time  at  which  the  blow  is  initiated.  C2  can 
vary  from  —0.1  to  —0.01  s_1  depending  on  nozzle  diameter  and  reservoir  size.  For  a  particular 
boat,  reservoir  size  depends  on  which  and  how  many  reservoir  options  the  operators  choose  to 
use. 

Some  analyses  simply  assume  the  mass  flow  rate  from  the  reservoir  is  proportional  to  the 
mass  left  in  the  reservoir  [16].  This  is  equivalent  to  setting  k  =  1  in  (42)  which  results  in: 


Prify  Pro ® 


C2  (t  to) 


(46) 


This  ‘exponential’  model  is  compared  to  the  isentropic  model  in  Figure  7.  With  identical 
blowing  constants,  the  exponential  model  empties  the  reservoir  faster  then  the  ideal  frictionless 
isentropic  model,  which  doesn’t  make  sense.  However,  using  a  smaller  blowing  constant  with 
the  exponential  model  provides  a  time  response  similar  to  that  of  the  isentropic  model.  The 
model  to  use  depends  on  how  the  blowing  constant  is  obtained.  If  the  blowing  constant  is 
calculated  from  a  known  initial  mass  and  mass  flow  rate,  then  the  isentropic  model  should  be 
more  accurate.  If  C2  is  obtained  by  fitting  time  response  data  to  the  exponential  model,  then 
the  exponential  model  should  be  used.  Herein,  the  exponential  model  is  used  in  keeping  with 
the  precedent  set  by  Mackay  [16]. 

The  total  mass  of  air  ma  in  the  MBTs  must  equal  that  released  by  the  reservoir  and  the 
sum  of  the  air  masses  mai  in  each  individual  tank: 

N 

ma  =  mro  (l  -  ^  mai  (47) 

i=  1 


It  is  now  assumed  that  the  air  delivery  system  to  the  MBTs  is  tuned  so  that: 


rngVri 

VT 


(48) 


The  next  assumption  is  questionable  and  is  made  to  avoid  having  to  use  a  heat  transfer 
model.  It  is  known  that  the  air  jetting  into  the  tanks  generates  good  mixing,  creating  a 
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Figure  7  Two  blowing  models  giving  similar  blow  histories  using  different  blowing  constants. 


turbulent  and  frothy  air /water  interface.  This  promotes  good  heat  transfer  from  the  water 
to  the  expanding  cold  air.  Therefore,  the  assumption  made  is  that  the  air  blowing  into  the 
MBTs  takes  on  the  ambient  MBT  temperature  immediately.  This  generates  buoyancy  more 
quickly  than  will  actually  occur  and,  in  the  current  context,  may  bring  on  the  rising  instability 
prematurely. 

The  air  pressure  in  ballast  tank  i  is  therefore: 

Pai  =  Pat  +  P9zwi  =  Pai.RT  (49) 


where  pat  is  atmospheric  pressure,  zwi  is  the  depth  of  the  water  level  in  tank  i  below  the 
ocean  surface,  pai  is  the  density  of  the  air  in  tank  i,  R  is  the  gas  constant  for  air,  and  T  is  the 
ambient  temperature  in  the  MBTs.  (49)  ignores  the  pressure  loss  through  the  openings  in  the 
bottom  of  the  MBTs  through  which  flood  water  is  expelled  during  the  blow.  These  openings 
are  large  so  the  MBTs  fill  quickly  when  the  boat  dives;  losses  for  the  slower  rising  maneuver 
should  be  small. 


The  volume  of  air  in  each  MBT  is  then: 


T/  _mai  _mro{l-ec*-^)RTVTi 

*  n.n 


(50) 


Pai  (Pat  +  PgZwi)VT 

When  the  boat  is  pitched  at  an  angle  6,  zwi  is  different  for  each  MBT  and  is  itself  dependent 


on  Vai: 


zwi  =  zo  ~  xTi  si11  $  —  0.45d  cos  0  1  — 


2K 

VT, 


(51) 


This  is  only  valid  for  zwi  >  0.  Substituting  (51)  into  (50)  results  in  an  expression  quadratic  in 
VajVTii  the  solution  for  which  is: 


^  =  A,  +  X/Af  +  A2 

VT  i 


(52) 


Ax  — 


—Pat  “  P9  (zo  “  xTi  sin  6  “  0.45d  cos  9) 
l.&pgdcos  9 


Ao  — 


m. 


■  0RT(1  _eCa(t-io)) 


0.9 pgdVT  cos  9 
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(a)  (b) 


Figure  8  Blowing  ballast  (a)  at  a  pitch  angle  of  20  degrees  at  a  fixed  depth  of  50  m  and  (b )  while 
rising  to  the  surface  from  an  initial  depth  of  100  m.  C2  =  —0.045  s~L,  mro  =  1500  kg,  the 
xTi  and  other  parameters  are  given  in  Appendix  A. 


with  the  constraint  that  Vai/VTi  is  limited  to  1.0.  When  Vai/VTi  >  1,  air  is  expanding  out 
the  bottom  of  the  MBT.  During  a  simulation,  zwi  should  be  checked  to  ensure  it  is  positive;  a 
negative  zwi  value  indicates  the  boat  is  on  the  surface. 

Figure  8  shows  the  predictions  of  the  blowing  model  for  a  boat  with  a  fixed  pitch  angle 
sitting  at  a  fixed  depth  in  one  case  (unrealistic)  and  coming  to  the  surface  at  a  constant  rate 
in  another.  In  the  latter  case,  the  boat  is  still  35  m  deep  when  air  begins  escaping  from  the 
forward  ballast  tank.  This  figure  shows  the  importance  of  modelling  depth,  pitch  angle,  and 
ballast  tank  axial  location. 

Finally,  it  is  clear  that  ballast  blowing  will  result  in  discontinuities  in  the  equations  of 
motion  at  the  times  that  the  MBTs  empty.  These  discontinuities  can  compromise  the  efficiency 
and  accuracy  of  the  ODE  integrator  if  they  are  severe  enough.  Predicting  the  times  at  which 
discontinuities  occur  so  the  integration  can  be  stopped  and  restarted  at  these  times  requires, 
from  (52),  knowledge  of  0  and  z0  at  the  discontinuities,  states  which  are  themselves  solutions 
of  the  integration  and  therefore  unknown  until  the  discontinuity  has  been  reached.  Methods 
are  available  for  handling  this  problem  but  are  not  implemented  at  this  time.  In  practice, 
integrating  through  the  discontinuities  has  not  proven  to  be  a  problem;  some  inaccuracy  results 
but  not  as  much  as  is  present  in  the  modelling  assumptions  themselves. 
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6  Control  System  Modelling 

When  control  surface  deflections  or  propulsion  changes  are  commanded,  the  changes  are  imple¬ 
mented  through  control  systems  of  varying  complexity.  Following  Watt  [17],  these  changes  are 
modelled  using  a  second  order  differential  equation: 

+  ^28  =  u2dc  (53) 

where  6  is  the  time  dependent  quantity  being  modelled  (control  surface  deflection  or  com¬ 
manded  forward  speed),  Sc  is  the  commanded  value  of  5  (ie,  the  value  the  control  system  is 
trying  to  implement),  (  is  control  system  damping  (assumed  sub-critical:  0  <  (  <  1),  and  uj 
is  the  response  frequency  of  the  control  system.  The  general  solution  to  (53)  is: 

8  =  8C-  ^gin^°e~Cu,(t~to)  sin  (V1  -  C 2v(t  -  t0)  +  (3s)  (54) 

where  t0  is  the  time  at  which  the  new  command  5C  is  issued,  50  is  the  value  of  5  at  t  =  t0,  and 
/3  is  a  phase  shift  used  to  match  the  rate  d5/dt  at  t  =  t0.  Thus,  S0  and  /3  ensure  the  model 
for  6  is  first  order  continuous  when  a  new  command  is  issued.  The  enabling  mathematics  and 
some  implementation  examples  are  provided  in  [17]. 

To  use  this  model,  three  invariant  control  system  characteristics  must  be  specified  for  each 
control  system: 

£  Control  system  damping  as  described  above.  The  lower  the  damping  the  faster 

the  system  achieves  5C,  but  at  the  expense  of  overshoot.  Typical  values  range 
from  0.7  to  0.9. 

tUmax  The  natural  response  frequency,  the  maximum  frequency  at  which  the  system 
can  respond. 

The  rate  limit,  the  maximum  rate  at  which  the  system  can  respond. 

RL 

The  response  to  any  given  command  is  determined  by  first  assuming  uj  =  cumax  and  then 
checking  to  see  if  the  response  rate  dS/dt  is  less  than  the  rate  limit.  If  true,  this  ‘frequency 
limited’  solution  is  implemented.  If  false,  a  ‘rate  limited’  solution  is  used  in  which  c o  is  reduced 
until  the  maximum  d5/dt  magnitude  in  the  response  matches  the  rate  limit. 

This  algorithm  can  be  implemented  using  the  two  FORTRAN  subroutines  listed  in  Ap¬ 
pendix  B.  The  first,  CNSYS2,  calculates  the  50,t0,/3,uj  parameters  at  the  issuance  of  each  new 
5C  command  for  each  control  system.  These  are  saved  and  passed  to  the  second  subroutine, 
CSDEFL,  which  simply  calculates  (54)  at  any  point  in  time. 

Figure  9  shows  the  control  system  model  in  action.  Except  for  different  damping  param¬ 
eters,  the  invariant  control  system  characteristics  are  the  same  in  each  half  of  the  figure,  as  is 
the  initial  rate  d5/dt  =  0  at  t  =  t0.  In  part  (a),  the  response  to  5C  =  5  is  frequency  limited, 
with  uj  =  u;max  whereas  the  response  to  6C  =  10  is  rate  limited  with  uj  =  1.25  s^1. 

In  Figure  9b,  new  commands  are  issued  every  second  for  the  first  three  seconds  to  show 
how  the  algorithm  maintains  first  order  continuity  in  its  response. 


d5 

dt 
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a)  Initial  commands  5C  =  5,10. 


b)  Changing  commands  at  t  =  0, 1, 2,  3  s. 


Figure  9  First  order  continuous  control  system  time  responses  to  different  initial  commands 
and  a  sequentially  changing  command.  tomax  =  2  s-1,  (dd/dt)^  =0.1  6  units  per  second. 


7  Estimating  Roll  Stability 

Watt  [4]  presents  a  1  DOF,  quasi-steady,  linearized  stability  analysis  of  (f>  in  the  rolling  moment 
equation  of  motion  assuming: 

•  9,q,r  =  0  and  p  =  (j), 

•  (f>  is  decoupled  from  0,Q,u  so  these  latter  parameters  that  can  be  arbitrarily  large, 

•  BG  is  constant. 

The  result  is  a  simple  expression  giving  the  balance  between  the  stabilizing  static  and  destabi¬ 
lizing  hydrodynamic  forces: 


BBG  cos  9 


\p^u2 


OK' 

~d¥ 


>  0 

<E>=0 


for  stability. 


(55) 


The  second  term  is  the  destabilizing  force  and  it  varies  with  the  dynamic  pressure  of  the  flow 
and  so  increases  with  the  square  of  the  velocity.  It  is  also  proportional  to  a  stability  derivative 
that  can  be  estimated  from  (A5): 


=  0.200  cos  0  sin  0  +  0.891  sin2  0  +  0.449  sin3  0  (56) 

$=o 

This  is  nonlinear  in  flow  incidence  0,  as  shown  in  Figure  10. 

Watt  [4]  discusses  how  the  quasi-steady  assumption  leads  to  an  underprediction  of  the 
stability  derivative  when  this  derivative  is  obtained  from  a  moment  measurement  in  a  steady 
flow  with  the  trailing  vortex  field  from  the  sail  fully  developed.  This  trailing  wake  interacts 
with  the  tailplanes  reducing  the  rolling  moment  generated  by  the  sail  alone.  The  problem 
with  the  steady  flow  data  is  that  when  a  boat  starts  to  roll,  the  trailing  vortex  field  is  not 
fully  developed;  indeed,  it  is  just  beginning  to  develop.  Hence,  there  is  some  justification  for 
using  stability  derivatives  from  steady  state  measurements  with  the  tail  removed.  Reference  [4] 
compares  the  two  methods.  Herein,  (56)  is  used  for  simplicity. 
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Figure  10  The  nonlinear  hydrodynamic  stability  derivative. 


Watt  and  Bohlmann  [5]  also  show  how  q  and  q  effect  the  roll  stability  limit.  For  simplicity 
again,  this  effect  is  neglected  here. 

In  the  simulations  presented  in  the  next  section,  (55)  provides  a  simple  way  to  estimate 
stability.  Stability  is  monitored  using  a  velocity  stability  ‘index’: 


Uo  = 


N 


BBG  cosO 


\pf 


die 


d<h 


-u  >  o 


for  stability 


<3>=0 


(57) 


This  gives  physical  meaning  to  nonzero  Ug  values:  when  Us  =  1  m/s,  say,  then  the  boat  would 
be  unstable  if  U  were  1  m/s  larger,  everything  else  being  the  same.  Note  that  BG  has  been 
replaced  with  BG  from  (40)  so  the  stability  index  accounts  for  the  temporary  increase  in 
static  stability  while  the  tanks  are  being  blown. 

It  is  worth  examining  the  impact  of  two  conventional  remedies  for  the  rising  instability 
on  the  stability  index.  The  first  remedy  suggests  that  increasing  speed  is  beneficial  because 
it  reduces  flow  incidence  0  (see  Figure  1).  This  may  have  a  net  benefit,  even  though  the 
destabilizing  hydrodynamic  force  increases  as  the  square  of  the  speed,  since  from  (56)  when  0 
is  large,  dK/d<h  can  decrease  as  the  cube  of  0  ~  —w/U.  Thus,  in  (55)  the  highest  order  term 
in  U2(dK/d&)  is  0(1/17).  This  remedy’s  dependence  on  a  third  order  term  in  the  stability 
derivative  suggests  that  any  benefit  is  likely  to  be  both  marginal  and  geometry  dependent. 

The  second  remedy  suggests  that  increasing  pitch  angle  is  beneficial  because  it  reduces  0  by 
shifting  the  impact  of  buoyancy  away  from  increasing  —w  towards  increasing  u  and  therefore 
speed.  If  —w  falls  off  as  cos0,  then  the  highest  order  term  in  U2  (dK/d&)  is  0(cos3  6 /U) 
and  the  second  order  term  would  also  decrease  as  O(cos20).  This  effect  is  moderated  by  the 
reduction  of  static  stability  with  cos  9  in  (55),  so  that  the  net  effect  on  the  second  and  third 
terms  is  0(cos9)  and  0( cos2  9/U)  which  is  still  better  than  that  provided  by  the  first  remedy. 
Because  this  benefit  is  seen  at  lower  order,  it  is  more  likely  to  be  generally  useful. 
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8  Simulations  Using  the  Coefficient  Based  Model 

Rising  simulations  are  now  presented  using  the  coefficient  based  hydrodynamic  model  in  place 
of  the  RANS  model.  Various  scenarios  are  examined  to  see  how  they  effect  the  roll  angle  of 
the  boat  immediately  before  surfacing  (the  emergence  roll  angle).  A  typical  scenario  is  given 
in  Figure  1  from  Watt  and  Bohlmann  [5]  and  is  reproduced  here  as  Figure  11. 

The  boat  in  Figure  11  is  at  depth  when  it  uses  sternplane  control  to  pitch  up  to  the  desired 
angle  6  for  the  rising  maneuver.  It  then  blows  only  half  of  its  MBTs.  The  downward  force 
from  the  sternplanes  and  the  forward  momentum  of  the  boat  initially  generate  a  positive  w 
velocity  before  the  blow  progresses  enough  that  buoyancy  makes  w  negative;  this  accounts  for 
the  early  jump  of  n  in  4>,  shown  in  the  figure  as  a  change  of  sign  in  0*.  Although  propeller  rpm 
is  constant  throughout  the  maneuver,  the  axial  velocity  u  increases  because  of  the  buoyancy 
component  in  the  axial  direction,  and  from  the  thrust  resulting  from  the  hull  ‘sailing’  in  the 
crossflow.  Flow  incidence  0  increases  as  the  blowing  progresses  but  tapers  off  after  blowing 
stops.  The  air  continues  to  expand  in  the  ballast  tanks,  but  incidence  is  kept  in  check  by  the 
increasing  forward  speed.  The  maneuver  was  nominally  carried  out  in  a  vertical  plane  (the 
rudder  was  fixed)  but  the  heading  ip  still  varies. 

Roll  angle  (j)  is  small  through  most  of  the  rise,  until  just  after  the  MBTs  empty  and 
expanding  air  in  the  tanks  has  begun  escaping  and  possibly  interacting  with  the  sail.  The 
boat  begins  to  roll  and  emerges  through  the  surface  with  a  small  to  moderate  roll  angle.  As 
it  surfaces,  the  submarine  immediately  loses  static  stability  which  is  regained  gradually  as 
floodwater  drains  from  the  sail  and  deck  casing.  This  loss  of  stability,  combined  with  the 
emergence  roll  angle  and  roll  rate,  result  in  an  excessive  roll  angle  before  the  boat  recovers. 

Despite  the  care  taken  conducting  this  maneuver  (pitch  up  attitude  prior  to  blowing,  mod¬ 
erate  blowing),  an  uncomfortably  large  roll  angle  still  occurred  on  the  surface.  The  temporary 
loss  of  hydrostatic  stability  a  submarine  experiences  when  surfacing  is  well  understood  and 
normally  not  a  problem;  however,  this  surface  instability  is  aggravated  by  emergence  roll  which 
is  the  result  of  the  underwater  roll  instability  that  previous  computer  simulations  have  not  been 
able  to  predict. 

With  this  scenario  in  mind,  nine  simulations  with  the  present  model  have  been  carried 
out  and  are  presented  in  detail  in  Appendix  C.  The  differences  between  the  simulations  are 
summarized  in  Table  1,  significant  values  for  some  states  are  listed  in  Table  2,  and  the  roll 
histories  are  all  compared  in  Figure  12. 

The  equations  of  motion  are  formulated  and  solved  using  Maple.1  Within  Maple,  the  ODEs 
are  numerically  integrated  using  RKF45,  a  Runge-Kutta  integrator  that  chooses  its  own  step 
size  At  based  on  local  error,  the  difference  between  the  4th  and  5th  order  accurate  solutions 
it  propagates  simultaneously.  Time  steps  are  about  half  a  second  for  most  of  the  integration 
but  are  typically  scaled  back  to  half  that  to  get  past  discontinuities  associated  with  the  ballast 
tanks  emptying. 

Simulation  1 

This  first  simulation  provides  a  baseline  rising  maneuver  with  which  subsequent  simulations 
are  compared.  The  procedure  is  typical  of  that  used  by  some  of  the  larger  diesel  boats.  Unlike 
the  procedure  used  in  Figure  11,  Simulation  1  (SI)  increases  commanded  speed  when  ballast 
is  blown,  a  common  practice  intended  to  both  reduce  flow  incidence  and  get  the  boat  to  the 
surface  quickly. 
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Figure  11  Full  scale  rising  trial  for  a  2000-3000  t  boat  blowing  limited  ballast  at  depth.  0*  is 
0  allowed  to  go  negative  (for  display  purposes  only)  when  $  jumps  by  it. 


SI  is  displayed  in  Figure  Cl  in  Appendix  C  (the  first  page  of  this  appendix  explains  how 
to  read  the  figure).  It  begins  with  the  boat  in  equilibrium  in  straight  and  level  flight  with  a 
speed  of  3  m/s  (just  under  6  knots)  at  a  depth  of  100  m.  At  t  =  0,  the  sternplanes  are  activated 
to  pitch  the  nose  up,  ballast  is  blown  using  normal  air  only,  and  a  speed  of  6  m/s  is  requested. 
As  the  large  forward  MBTs  empty,  the  sternplane  deflection  is  reduced  to  keep  the  pitch  angle 
at  about  20  degrees. 

Figure  Cla  shows  that  u  initially  decreases.  This  is  because  of  increased  drag  from  the 
deflecting  sternplanes.  After  10  to  15  seconds  u  increases  quickly  as  uc  and  the  pitch  angle  and 
net  buoyancy  increase.  At  the  end  of  the  maneuver,  u  is  still  increasing  because  of  buoyancy 
and  the  nose  up  pitch  angle.  The  crossflow  —w  increases  steadily  as  the  buoyancy  increases. 

The  pitch  rate  q  (Figure  Clb)  initially  increases  because  of  the  sternplane  deflection  but 
levels  off  as  sternplane  control  is  countered  by  static  stability  in  pitch.  It  suddenly  decreases 
when  the  sternplane  deflection  is  reduced  to  maintain  the  pitch  angle,  q  eventually  increases 
again,  forcing  the  sternplanes  to  stay  active,  because  of  increasing  buoyancy  in  the  nose  MBTs 
which,  as  well  as  being  larger,  are  under  less  pressure  than  the  stern  tanks.  This  stops  suddenly 
when  the  forward  tanks  empty,  after  which  increasing  buoyancy  in  the  aft  MBTs  reduces  q. 

As  the  commanded  speed  uc  in  Figure  Cla  increases  from  3  to  6  m/s,  p  in  Figure  Clb  and 
4>  in  Cld  increase  in  magnitude  as  propeller  torque  increases  (Jb  decreases).  4>  is  more  clearly 
seen  in  Figure  12.  Since  buoyancy  does  more  to  accelerate  the  boat  than  propulsion,  Jb  begins 
to  increase  at  about  t  =  18  s  and  this  will  cause  propeller  torque  to  decrease  as  well.  Despite 
decreasing  torque,  both  p  and  4>  steadily  increase  in  magnitude  for  the  rest  of  the  maneuver  as 
the  boat  responds  to  the  increasing  moment  on  the  sail. 

The  mass  characteristics  are  shown  in  Figure  Cle.  The  discontinuity  that  results  from  the 
forward  ballast  tanks  emptying  generates  discontinuities  in  all  the  mass  characteristics  but  does 
not  cause  the  ODE  integrator  much  trouble,  probably  because  the  effect  is  small  relative  to 
the  overall  mass  of  the  vehicle.  As  discussed  earlier,  the  static  stability  BG  increases  steadily 
through  the  maneuver  and  then  falls  relatively  quickly  at  the  end,  just  as  the  tanks  empty,  to 
its  original  value.  The  discontinuity  as  the  tanks  empty  is  an  artifact  of  the  current  model  and 
could  be  eliminated  by  refining  the  model. 
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Simulation 

Initial  Conditions:  unless  otherwise  noted,  equilibrium  flight  at  uc  =  3  m/s,  <j>  =  <f>0, 

5S  =  5st,  zo  =  100  nr,  xG,yG  =  0;  additional  parameters  as  per  Appendix  A. 

6s  uc 

commanded  commanded 

(deg.)  (m/s)  Blow  Characteristics 

SI 

— 20  @  1  =  0  6@1  =  0  normal  Baseline  simulation:  4>0  =  —0.142  degrees, 

—  1  @1  =  17  (5S t  =  0.813  degrees,  yGo  =  0, 

3  @  1  =  38  Wt  =  lM(10~4)B. 

S2 

— 20@1  =  0  3  normal  SI  without  the  speed  increase. 

-1  @1  =  17 

3  @  t  =  38 

S3 

— 11  @  t  =  0  4.5  @  t  <  0  normal  SI  with  50%  more  speed,  twice  the  dynamic 

—  1  @1  =  17  9@1  =  0  pressure.  Initial  equilibrium  at  uc  =  4.5  m/s 

1  @  t  =  33  =>-  (j)0  =  0.320  deg.,  Wt  =  3.69(10“4)B. 

S4 

6S t  6@1  =  0  normal  SI  without  using  the  sternplanes. 

S5 

—20  @1  =  0  6  @  1  =  0  emergency  SI  except  the  MBTs  are  blown  using  the 

1  @  t  =  16.5  emergency  air  which  provides  twice  the 

3  @  t  =  35  mass  of  air  as  for  a  normal  blow. 

S6 

— 20  @  f  =  0  6@f  =  0  normal  SI  except  when  the  boat  reaches  a  depth  of 

—  1  @t  =  17  50  m,  the  sternplanes  are  used  to  reduce  the 

25  @  t  =  37.7  pitch  angle  to  about  5  degrees  on  surfacing. 

S7 

—25 @t  =  0  6@t  =  0  normal  SI  except  the  boat  is  surfaced  at  a  higher 

1  @  t  =  25  pitch  angle  of  35  instead  of  20  degrees. 

S8 

—20 @1  =  0  6@1  =  0  normal  SI  except  the  magnitude  of  the  initial  heel 

—  1  @1  =  17  angle  <p0  is  increased  to  2  degrees  by  setting 

3  @  1  =  38  yGo  =  -0.0114  m. 

S9 

—9 @1  =  0  6@1  =  0  normal  S8  except  the  sternplanes  limit  the  pitch 

0  @  1  =  16  angle  to  about  10  degrees  throughout 

3  @  1  =  40  the  maneuver. 

Table  1  Eight  variations  on  a  baseline  simulation  (SI ).  The  simulations  are  shown  in  Appen¬ 
dix  C  and  use  the  coefficient  based  hydrodynamic  model  from  Appendix  A. 

The  normal  velocity  w,  while  it  is  never  positive,  has  a  maximum  at  about  t  =  4  s  because 
of  the  downward  force  from  the  sternplanes  and  forward  momentum  of  the  boat.  Because  v  is 
also  small,  $  is  very  sensitive  to  this  (Figure  Clf)  but  eventually  settles  down  as  w  becomes 
decisively  negative. 

The  stability  index  Us  is  initially  high  and  positive  but  drops  rapidly  as  speed  and  flow 
incidence  increase.  Stability  is  lost  just  before  the  forward  tanks  empty,  when  BG  is  decreas¬ 
ing.  This  timing  is  consistent  with  that  of  the  trial  in  Figure  11.  Despite  the  instability,  roll 
magnitudes  are  not  large  enough  to  be  of  concern. 
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Simulation 

d7>* 

^max 

^  *  b 

Icq 

=  0 

tu 

te 

te 

(s) 

Emergence  States 

Roll 

Angle 

Us 

(m/s)  (deg.) 

(at  t  — 

Pitch 

Angle 

6 

(deg.) 

te) 

Flow 

Incidence 

0 

(deg.) 

BGo 

BG0 

SI 

1.213 

1.15 

0.92 

47.4 

-1.81 

-2.0 

20.5 

17.1 

S2 

1.213 

1.17 

0.91 

48.6 

-1.75 

0.3 

21.5 

19.4 

S3 

1.213 

1.08 

0.94 

43.0 

-1.46 

-3.6 

20.4 

12.0 

S4 

1.217 

1.21 

0.86 

62.1 

-2.24 

-3.9 

9.5 

21.5 

S5 

1.215 

1.16 

0.83 

43.6 

-2.68 

-4.1 

20.1 

17.7 

S6 

1.213 

1.19 

0.86 

49.2 

-2.58 

-3.4 

5.4 

29.7 

S7 

1.209 

1.08 

0.97 

42.5 

-0.55 

-1.3 

35.8 

12.1 

S8 

1.213 

1.15 

0.91 

47.5 

-1.82 

-9.2 

20.1 

17.2 

S9 

1.216 

1.20 

0.86 

55.0 

-2.26 

-16.4 

10.2 

21.1 

Table  2  Some  simulation  results:  The  maximum  BG  value  is  virtually  constant.  BGmax  and 
BGu,  the  BG  value  just  as  the  boat  becomes  unstable,  are  shown  relative  to  BG0  =  0.35  m. 
The  time  at  which  the  boat  becomes  unstable,  tjj,  is  normalized  by  the  time  taken  to  surface. 
The  remaining  values  are  those  at  the  end  of  the  maneuver,  at  t  =  te. 
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Simulation  2 

52  differs  from  SI  only  by  maintaining  a  constant  commanded  speed  of  3  m/s  throughout  the 
maneuver.  This  reduces  u  by  only  about  10%  at  the  end  of  the  maneuver  since  u  is  strongly 
influenced  by  buoyancy  and  the  pitch  angle.  Flow  incidence  at  the  surface  is  increased  from 
17  to  19  degrees  (Table  2).  The  net  effect  is  slightly  detrimental.  There  is  a  slight  decrease  in 
stability:  BGjj  is  slightly  larger  so  instability  has  occurred  sooner  in  the  blow  (tjj  relative  to 
the  time  at  which  the  Val/Vn  reaches  1.0).  As  well,  there  is  slightly  more  time  for  the  boat 
to  respond  to  the  instability  (trj /tc  is  slightly  smaller). 

The  most  interesting  difference  between  SI  and  S2  is  the  effect  on  p  and  <j).  Since  uc  is 
constant,  the  advance  ratio  initially  increases  (it  initially  decreases  in  SI)  and  eventually 
increases  by  a  factor  of  three  as  u  increases.  This  results  in  decreasing  torque  and  a  decreasing 
heel  angle;  p  magnitudes  are  always  positive  and  so  the  boat  rolls  in  the  opposite  direction 
than  in  SI.  Thus,  while  propeller  thrust  has  only  a  minor  effect  on  this  maneuver,  torque  can 
influence  the  direction  of  roll.  One  wonders  if  the  propeller  could  be  used  to  maintain  a  constant 
roll  angle. 

Again,  the  roll  magnitudes  here  are  not  significant. 

Simulation  3 

53  differs  from  SI  by  increasing  the  commanded  speeds  by  50%  which  doubles  the  dynamic 
pressure  (hydrodynamic  forces)  throughout  the  maneuver.  Sternplane  control  is  adjusted  to 
keep  about  the  same  maximum  pitch  angle  of  20  degrees.  Despite  the  faster  speed,  the  boat 
reaches  the  surface  only  4.4  s  sooner  than  in  SI.  This  is  because  little  depth  change  occurs 
in  the  first  half  of  either  maneuver  and  —  w,  about  half  of  the  rising  velocity  (cf,  Figure  1),  is 
slightly  reduced  in  S3  relative  to  SI. 

A  substantial  reduction  in  0  occurs  over  most  of  this  maneuver  because  of  the  large  increase 
in  u  and  slight  decrease  in  —w.  Although  this  is  countered  to  a  large  extent  by  the  increase 
in  dynamic  pressure,  BGjj  shows  that  instability  has  been  delayed.  This  is  also  seen  in  the 
increased  tu/te  ratio. 

S3  brings  to  light  a  disturbing  new  parameter  that  needs  to  be  considered.  Despite  the  fact 
that  S3  is  the  shorter,  more  stable  maneuver,  it  achieves  a  roll  angle  80%  greater  than  in  SI. 
This  is  attributed  to  the  larger  initial  heel  angle  (c f>0  =  —0.32  versus  —0.14  degrees  in  SI)  which 
results  from  the  higher  initial  propulsion  torque  required  for  the  S3  maneuver.  Clearly,  initial 
heel  is  an  important  factor  in  this  analysis. 

The  S3  emergence  roll  angle  is,  nevertheless,  not  large  enough  to  be  of  concern. 

Simulation  4 

54  differs  from  SI  in  that  sternplane  control  is  not  used  to  pitch  the  nose  up.  This  results 
in  a  much  longer  maneuver.  The  crossflow  —w  is  primarily  determined  by  the  buoyancy;  it 
develops  earlier  in  the  absence  of  downwards  force  from  the  sternplanes  but  not  as  quickly  as  in 
SI  because  the  ballast  tanks  empty  slower.  However,  —w  achieves  a  higher  magnitude  because 
it  has  more  time  to  develop.  This,  combined  with  a  lower  axial  velocity,  results  in  substantially 
larger  flow  incidence  throughout  the  maneuver. 

The  net  effect  is  a  measurable  decrease  in  stability,  as  seen  by  the  larger  BGjj  value  and 
substantially  lower  tjj/te  ratio.  Figure  12  shows  that  the  lower  pitch  angle  in  the  S4  maneuver 
has  little  effect  on  the  roll  time  history  prior  to  the  SI  maneuver  becoming  unstable.  The  big 
difference  between  the  simulations  is  the  extra  time  the  S4  maneuver  gives  the  boat  to  respond 
to  the  instability  once  it  occurs.  This  results  in  a  roll  angle  almost  twice  that  of  SI. 
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The  S4  emergence  roll  angle  is  still  not  large.  The  Figure  11  emergence  roll  angle  is  much 
larger  without  the  aggravating  factor  of  an  extremely  low  pitch  angle. 

Simulation  5 

55  is  the  same  as  SI  except  that  emergency  ballast  blowing  is  used  instead  of  normal  blowing. 
This  means  that  twice  the  mass  of  air  is  blown  in  S5  than  in  SI,  resulting  in  both  fore  and 
aft  tanks  emptying  before  the  boat  surfaces.  Slightly  different  sternplane  control  is  used  to 
maintain  a  pitch  angle  of  20  degrees.  The  maneuver  occurs  about  4  s  faster  than  SI  because 
full  buoyancy  is  achieved  faster.  Although  —w  grows  faster  than  in  SI,  it  is  largely  countered 
by  higher  u  magnitudes  so  0  is  only  slightly  larger.  BGjj  is  little  changed  but  tv/te  is  the 
smallest  of  all  the  maneuvers,  even  though  te  is  also  small.  Thus,  the  extra  buoyancy  gets  the 
boat  to  the  surface  in  92%  of  the  time  of  the  SI  maneuver  but  also  gives  the  boat  twice  as  much 
time  to  react  to  the  instability,  resulting  in  twice  the  roll  angle. 

Simulation  6 

56  is  the  same  as  SI  except  the  sternplanes  are  used  to  curtail  the  pitch  angle  as  the  boat 
passes  the  50  m  depth  mark  so  the  boat  can  surface  on  a  relatively  even  keel.  This  generates 
higher  flow  incidence  at  the  end  of  the  maneuver  which  is  partially  countered  by  lower  dynamic 
pressure.  It  brings  on  the  instability  somewhat  earlier  and  delays  completion  of  the  maneuver, 
increasing  the  time  during  which  the  boat  is  unstable  by  80%.  This  results  in  a  70%  increase 
in  the  final  roll  angle  although,  again,  the  final  roll  angle  is  not  large. 

Full  scale  maneuvers  suggest  that  last  minute  pitch  curtailment  can  result  in  substantial 

roll. 

Simulation  7 

57  is  the  same  as  SI  except  the  sternplanes  are  kept  deflected  longer  so  the  boat  surfaces  at 
a  pitch  angle  of  35  degrees.  This  action  is  opposite  in  nature  to  the  pitch  curtailment  used 
in  S6.  This  is  the  fastest  maneuver,  even  beating  S3.  Flow  incidence  is  substantially  reduced 
and  instability  is  much  delayed.  The  final  roll  angle  is  substantially  less  than  for  SI. 

Simulation  8 

58  is  the  same  as  SI  except  that  the  CG  is  moved  laterally  off  the  centerline  by  11  cm  generating 
an  initially  large  heel  (4>0)  of  2  degrees.  This  results  in  little  or  no  change  to  most  time  histories, 
including  flow  incidence,  speed,  and  the  stability  index.  However,  flow  incidence  now  generates 
a  much  larger  force  on  the  rolled  sail  (cf,  Figure  10)  which  produces  higher  roll  rates  and  an 
emergence  roll  angle  of  9  degrees. 

Both  the  initial  heel  and  final  roll  angles  here  are  similar  to  those  in  Figure  11.  The  final 
roll  clearly  occurs  as  a  result  of  the  relatively  large  initial  heel  angle. 

Simulation  9 

59  is  S8  with  the  pitch  angle  limited  to  about  10  degrees  throughout  the  maneuver.  This  brings 
on  instability  sooner  (larger  BGjj  value)  and  the  time  during  which  the  boat  is  unstable  (1 

is  doubled,  resulting  in  an  80%  larger  emergence  roll  angle.  This  is  very  much  a  concern  given 
the  consequences  upon  surfacing  shown  in  Figure  11. 
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Discussion 

Simulations  1  through  7  result  in  roll  angles  that  are  small  enough  to  be  significantly  influenced 
by  propeller  torque  changes.  This  complicates  the  interpretation  of  the  simulations.  Since  the 
roll  angles  are  too  small  to  be  of  concern,  the  conclusions  drawn  from  these  simulations  are 
tentative;  they  are: 

•  increasing  propeller  speed  has  a  small  beneficial  effect  on  roll  stability  and  is  important  for 
minimizing  the  duration  of  the  rise. 

•  propeller  torque  has  the  potential  for  growing  or  reducing  the  initial  heel  angle. 

•  buoyancy  is  necessary  but  excessive  buoyancy  unnecessarily  aggravates  the  roll  instability. 

•  increasing  pitch  is  the  best  way  to  increase  both  speed  and  roll  stability. 

•  pitch  curtailment  when  surfacing  aggravates  the  roll  instability. 

58  shows  that  a  large  (but  still  realistic)  initial  heel  angle  has  the  greatest  effect  on  the 
emergence  roll  angle.  It  does  not  effect  stability  because  static  stability  and  the  destabilizing 
hydrodynamic  force  both  increase  linearly  with  </>.  However,  the  rolling  moment  induced  by 
the  crossflow  grows  with  the  heel  angle,  and  the  additional  heel  needed  to  offset  this  additional 
moment  grows  as  well.  This  is  easily  seen  from  a  simplified  quasi-steady  analysis  of  the  rolling 
moment  equation.  If  is  small  so  $  ~  </>,  Q  is  propeller  torque,  and  yB  =  zG  =  0,  then  the  K 
equation  can  be  written  as: 

small  terms  =  K$<fi  —  Q  +  yGW  +  zBB(j) 

so  that: 

Q  —  yGW  +  (small  terms) 

Ii,j>  +  zbB 

Thus,  heel  grows  as  the  magnitude  of  the  denominator  in  (58b)  decreases;  that  is,  as  stability 
decreases.  This  is  why  </>  grows  in  the  simulations  even  though  the  boat  is  still  stable  —  new 
equilibria  are  being  sought  to  offset  the  increasing  rolling  moment.  However,  the  growth  in  cj) 
is  also  proportional  to  the  numerator  which  increases  with  the  initial  heel  angle  (determined 
by  Q  —  yGIH).  So  S3  and  S8  see  greater  roll  changes  prior  to  instability  than  does  SI.  When 
the  boat  goes  unstable  in  S3  and  S8,  it  does  so  with  a  higher  roll  rate  and  rolling  moment  so 
subsequent  roll  is  more  severe  than  for  SI.  (Note  that  (58)  becomes  invalid  when  instability 
occurs  because  the  dynamics  are  no  longer  quasi-steady.) 

59  confirms  the  finding  with  the  lower  amplitude  roll  simulations  that  decreased  pitch 
does  not  substantially  alter  the  roll  time  history  prior  to  instability.  Decreased  pitch  decreases 
stability  and  provides  increased  time  for  the  boat  to  react  to  the  instability.  This  increases 
the  emergence  roll  angle  and  roll  rate  which,  in  turn,  aggravate  roll  as  the  boat  surfaces  and 
temporarily  loses  static  stability. 
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9  Concluding  Remarks 

Propulsion,  blowing,  and  control  system  models  have  been  presented  and  tested  prior  to  use  in 
an  unsteady  6  DOF  RANS  simulation.  The  models  were  tested  using  a  high  incidence,  quasi¬ 
steady,  coefficient  based  hydrodynamic  model  and  several  submarine  emergency  rising  scenarios 
were  explored. 

The  coefficient  based  simulations  show  that  the  boat  consistently  becomes  unstable  just 
before  the  ballast  tanks  empty,  several  seconds  before  surfacing.  The  best  way  of  dealing  with 
the  roll  instability  is  to  use  speed,  induced  by  moderate  buoyancy  and  a  substantial  pitch  angle 
so  instability  is  not  aggravated,  to  get  the  submarine  to  the  surface  before  it  has  time  to  respond 
to  the  instability. 

When  the  initial  heel  angle  is  small  (much  less  than  a  degree),  emergence  roll  angles  are 
also  small.  However,  when  the  initial  heel  angle  is  appreciable  (eg,  2  degrees),  emergence  roll 
angles  are  a  concern;  combined  with  low  pitch  angles,  emergence  roll  angles  large  enough  to 
instigate  excessive  roll  after  surfacing  can  result. 

Further  investigations  are  needed  to  validate  these  preliminary  findings.  The  unsteady 
RANS  simulations  to  follow  will  be  important  in  this  regard.  Full  scale  trials  are  necessary  to 
validate  the  suggested  mitigation  strategies  and  to  see  what  kind  of  heel  angles  boats  typically 
run  with.  The  isothermal  blowing  model,  which  assumes  the  blown  air  temperature  immediately 
warms  to  ambient  temperatures,  needs  further  investigation. 
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Appendix  A:  Simulation  Constants,  Functions,  and  Coefficients 


Constants 

The  following  constants  are  used  in  the  simulations: 


k2mps  =  0.514 
g  =  9.81 
t  =  70 
d  =  €/8.75 
D  =  d/2 
p  =  1025 
V  =  0.008645P3 


Convert  knots  to  m/s. 

Gravitational  constant,  m/s2. 

Overall  length  of  boat,  m. 

Maximum  hull  diameter,  m. 

Propeller  diameter,  m. 

Density  of  sea  water,  kg/rn3. 

Volume  displaced  by  the  hydrodynamic  envelope,  m3. 


BG0  =  0.005£ 


Initial  height  of  CB  above  CG,  m. 


xBi  Vbi  zb  —  0’  0,  BG0 

XGo >  VGoi  ZGo  =  XBiVBi  zb  +  BG0 

ho  =  3  (10-')  (\pf) 

4„  =  10“3  (PD 


ho- 

II 

O 

Ixy0 

’  dyz„ 

!  dZXQ  0. 

o 

cT 

N  = 

:  4 

Vti-, 

l/r2! 

Vt3 ,  Vt 4 

=  100,60,60, 

30 

VT  = 

=  Vti 

+  Ft2  + 

vT3  +  vT4  = 

250 

XT 1  • 

xT2i 

xT3,  xT4 

=  22,19,-29 

,-32 

T  =  278.16 
pa  =  101325 
R  =  287 

C2  =  -0.06,  -0.03 
mro  =  1000,  2000 

C,wmax,  dd/dt\RL  =  0.9, 2.0s-1, 0.1s_1 
C,^max,  dS/dt\RL  =  0.7,0. Is-1, 0.5m/s2 


CB  coordinates,  m. 

Initial  CG  coordinates,  m. 

Initial  moment  of  inertia  about  x  axis,  kg  m2. 
Initial  moment  of  inertia  about  y  axis,  kg  m2. 
Initial  moment  of  inertia  about  z  axis,  kg  m2. 
Initial  products  of  inertia,  kg  m2. 

Number  of  main  ballast  tanks. 

MBT  volumes,  m3. 

Total  volume  of  all  MBTs,  m3. 

Axial  location  of  MBT  centroids,  m. 

MBT  temperature,  degrees  Kelvin,  K. 
Atmospheric  pressure,  N/m2. 

Gas  constant  for  air,  J/kg/K. 

Normal,  emergency  blowing  constants,  s-1 
Normal,  emergency  air  mass  for  blowing,  kg. 
Sternplane  control  parameters. 

Propulsion  (uc)  control  parameters. 


Coefficient  Model  Translational  Hydrodynamic  Functions 

The  F,uvw  ( 0 ,  4>)  functions  used  in  the  coefficient  model  and  plotted  in  Figure  2  in  the  main  text 
are  listed  below. 


m)X'uvw  =  — 0.1460444814  cos2  (0)  +  2.092727148  sin2  (0) 

+  0.8822237467 sin4 (0)  +  0.1130677803  cos (0)  sin(0)  cos(<F) 

-  [0.7378986666 sin2 (0)  -  1.657751309 sin3(0)]  cos($) 

-  [0.8340902069 sin2 (0)  +  1.978042538 sin3(0)]  cos(2<4>) 

+  [0.04793423102 sin2(0)  -  0.3341608580 sin3 (0)]  cos(34>) 

-  [0.2580423873 sin2 (0)  -  1.195335953 sin3(0)]  cos(4<F)  (Al) 
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100 Y^vw  =  6.742784020  cos  (0)  sin(0)  sin(<4>) 

+  [11.01008861  sin2(0)  +  4.193261671  sin3(0)]  sin(4>) 

-  [12.89386476  sin2 (0)  -  20.17874874  sin3 (0)]  sin(2<4*) 

-  [0.3810901693 sin2(0)  -  2.091386770 sin3(0)]  sin(3$) 

-  [2.870012312 sin2(0)  -  8.508595821  sin3(0)]  sin(4$) 

-  [0.08969245281  sin2(0)  +  1.145O93633sin3(0)]  sin(54>)  (A2) 

100  Z'uvw  =  0.05545234430  cos2  (0)  +  5.718214158  sin2  (0) 

-  16.10350880 sin4(0)  +  1.984938259 cos (0)  sin(0)  cos($) 

+  [11.33288282 sin2(0)  -  33.70003698 sin3 (0)  +  47. 14241947 sin4 (©)]  cos($) 

-  [9.210462863 sin2(0)  -  15.58771719 sin3(0)  +  6.913780522 sin4(0)]  cos(2<4>) 

_  [4.447468749  sin2 (0)  -  30.70824244  sin3 (0)  +  36.56417687  sin4 (0)]  cos(3<4*) 

+  [3.004351162 sin2(0)  -  15.78083491  sin3 (0)  +  23.80409862 sin4 (0)]  cos(44>) 

+  [2.393554689 sin2 (0)  -  10.88133361  sin3(0)  +  9.034353557 sin4 (0)]  cos(5$)  (A3) 

lOOiv^^  =  0.2004601495  cos (0)  sin(0)  sin(4>) 

+  [0.9563341051  sin2(0)  -  0.7780718288 sin3(0)]  sin(<4>) 

+  [0.2673488329  sin2 (0)  -  0.1491263433  sin3  (0)]  sin(23>) 

-  [0.01868084765  sin2  (0)  -  0.1247368635  sin3 (0)]  sin(3$) 

-  [0.1360269435 sin2 (0)  -  0.2877553519 sin3(©)]  sin(44>)  (A4) 

100M',„,UJ  =  0.01548372486  cos2  (0)  +  1.024848874  sin2  (0) 

-  4.832975106 sin4(0)  -  0.4018487745  cos(0)  sin(0)  cos($) 

-  [0.1346352455 sin2 (0)  -  1.278018136sin3(©)]  cos($) 

-  [1.930064018 sin2(0)  -  3.352906727 sin3 (0)]  cos(2<4*) 

+  [0.2547988257 sin2 (0)  +  0.5809402243 sin3(0)]  cos(3$) 

+  [0.6306394660 sin2 (0)  -  0.5989226878 sin3(0)]  cos(44>) 

+  [0.4984216470  sin2 (0)  -  1.729231462  sin3 (©)]  cos(54>)  (A5) 

100KVW  =  1.409218893 cos(0)  sin(0)  sin(<3?) 

-  [2.020221148 sin2(0)  -  0.7333665071  sin3(0)]  sin(4>) 

+  [2.948685948  sin2 (0)  -  6.376794178  sin3(0)]  sin(2<4>) 

-  [0.1218281303 sin2(0)  -  0.2698840669 sin3(0)]  sin(3$) 

+  [0.3267021655 sin2(0)  -  1.560874539 sin3 (0)]  sin(4<4*) 

+  [0.1230668903 sin2 (0)  -  0.1030152504 sin3(0)]  sin(5$)  (A6) 
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DSSP20/DERIVS  Coefficient  Estimates 

The  DERIVS  program  [5]  uses  the  DSSP20  program  to  estimate  first  and  second  order  deriva¬ 
tives  for  a  simplified  submarine  geometry.  The  geometry  is  shown  in  Figure  Al.  The  control 
derivative  estimates  are: 


=  -0.0208151 
Z's  =  -0.0238569 

us 

Mi  =  -0.0109145 

us 

MSsSs  =  0 

mL\6s\  = 0 


X'5rSr  =  -0.0208151 
Y[r  =  0.0238569 
K'Sr  =  0 
^6,.  6,.  =  0 
N's  =  -0.0109145 


m 


DERIVS  generates  the  following  rotational  derivatives.  If  they  are  zero  or  very  small,  they  are 
replaced  with  their  added  mass  equivalent  if  that  is  appreciable.  The  numbers  shown  below 
are  the  dimensionless  coefficients  (eg,  X'pp  =  —0.0000281)  so  that  the  dimensions  of  the  (A8) 
coefficients  are  the  same  as  the  dimensions  of  the  terms  in  brackets  (eg,  the  dimensions  of  Xpp 
are  the  same  as  those  of  p£A). 


Xpp  =  -0.0000281  {^pl4) 

Xrr  =  0.0032739  (\pHA) 

Xwp  =  -0.0002454  (\pl3) 

Xuq  =  0.0000191  (ipf3) 

Xqq  =  0.0033202  (\ptA) 

Xq]ql  =  0.0000005  (\ptA) 

Yup  =  -0.0049919  (±p£3) 
YpM  =  -0.0000026  (\plA) 

Ywp  =  added  mass  equivalent 
Yur  =  0.0110609  (ipf3) 

Yr |r|  =  0.0030937  (^p£A) 

Zpp  =  added  mass  equivalent 
Zrr  =  0 

Zwp  =  -0.0000128  dpi3) 

Zuq  =  -0.0131470  (\pl3) 

Zqq  =  added  mass  equivalent 
Zq]q |  =  -0.0030973  (\pl4) 


Kup  =  -0.0004229  (ip£4) 
Kplp |  =  -0.0000002  {\p(.b) 
I\wp  =  added  mass  equivalent 
Kur  =  -0.0001373  (ip£4) 
Kr]r |  =  -0.0000002  (ip£5) 

Mpp  =  0.0000026  (ip£5) 

Mrr  =  0.0000031  (ip£5) 

Mwp  =  -0.0000033  (\ptA) 
Muq  =  -0.0060731  (ip£4) 

Mqq  =  0 

Mq]ql  =  -0.0012988  (ip£5) 

Nup  =  -0.0006396  {\p(.4) 
Np\p\  =  -0.0000003  (±p£5) 
Nwp  =  -0.0000248  (\pt4) 

Nur  =  -0.0064560  [\pt4) 
Nr\r\  =  -0.0012988  {\p(.b) 


(A8) 
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Figure  A1  The  simplified  geometry  used  by  DSSP20  [4]  for  estimating  vehicle  steady  state 
derivatives. 


Figure  A2  Planform  and  elevation  views  of  the  generic  shape  being  modelled  (black  lines)  and 
the  component  replacement  ellipsoids  (red  and  blue  lines)  used  for  estimating  vehicle  added 
mass,  following  [3]. 
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Added  Mass  Coefficients 

Replacing  vehicle  components  with  equivalent  ellipsoids  (Figure  A2),  for  which  analytical  ex¬ 
pressions  for  the  added  masses  are  available,  the  added  mass  coefficients  for  the  entire  vehicle 
can  be  estimated.  Interference  effects  on  the  appendages,  due  to  the  presence  of  the  hull  are 
estimated  [3].  The  coefficients  are  (cf,  Equations  (5)  &  (6)  from  the  main  text): 


Xi  =  -5.489  (10“4)  (\pl3) 

II 

X*  =  -2.844  (l0"6)(±p£3) 

Kp  =  -2.682  (10“5)  (ip£5) 

Xg  =  9.443  (10“6)  (\pt) 

Kf  =  -2.051  (lO"5)  (\plb) 

Yi}  =  -1.959  (KT2)  (lp£3) 

£ 

II 

Yp  =  -2.915  (KT4)  (\pt) 

.  -O' 

N 

II 

Yf.  =  2.067  (10“5)  (|p£4) 

Mg  =  -9.209  (10“4)  (ip£5) 

y  m  _  x  . 

II 

Zw  =  —1-622  (lO-2)  (|p^3) 

£ 

II 

=  -2.575  (KT4)  (\pt) 

Nf  =  -9.392  (10“4)  (ip£5) 

The  added  mass  coefficients  also  give  the  steady  state  potential  flow  forces  on  the  vehicle.  These 
are  only  used  when  there  is  no  estimate  available  that  accounts  for  viscous  effects.  The  steady 
state  potential  flow  coefficients  used  are  (see  Equation  25,  from  [3]): 


Y 

y*-vr 

X 

y*-rp 

X 

yvwq 


=  -Y,:, 
=  -Yp 
=  Zw 


V  —  x  ■ 

1  wr  ■/'-u 


y  —  v . 

i  qr  y\-q 

y  —  —  7 

1  wp  Z-Jyj 

Y  —  —7 

1  pq  zjq 


y 

z^wq 

Zqq 

y 

^vp 

y 

zjpp 

y 

£jrp 


=  -X; 


=  -x 

=  Vr 


K„ 


wp 


K 


pq 


K, 

K„ 


vq 


K, 


qr 


=  ~Y* 

=  Kf. 

=  Y,f.  +  Ztj 
=  ~Y,  -  Z* 
=  —Mg  +  Nj. 


MWq  =  Xg 

Mvr  =  Yp 
Mvp  =  -Y+ 

Mpr  =  Kp  -  Np 


N 

1  y  qr 

N 

1  y  vq 

N 
1  y  pq 


=  -K+ 

=  -X,  -  Yf  , 

=  ~Kp  +  Mg 


(A10) 
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Appendix  B:  Subroutines  for  Control  System  Modelling 

The  FORTRAN  77  program  cs.f  is  listed  below.  It  consists  of  a  main  program  and  two 
subrouines  CNSYS2  and  CSDEFL.  The  main  program  is  just  a  utility  for  exercising  the  sub¬ 
routines.  The  enabling  mathematics  for  the  subroutines  is  presented  in  [17]. 

PROGRAM  CS 


*  Front  end  for  SUBROUTINE  CNSYS2.  Reads  the  following  data  from  * 

*  input  file  SIMU.INP.  * 

*  * 

*  TEND  The  end  of  the  simulation.  * 

*  WMAX , DDMX , ZETA  Control  system  charateristics : response  * 

*  <leave  no  space>  frequency,  rate  limit,  damping.  * 

*  NCOM  Number  of  commands  to  be  issued  (max  20) .  * 

*  TC(1),DC(1)  Time  command  is  issued,  Defltn  commanded.  * 

*  TC (2) , DC (2)  ...  ...  * 

*  * 

*  TC (NCOM) , DC (NCOM)  * 

*  * 

*  Program  CS  calculates  the  proper  time  response  to  this  series  of  * 

*  commands,  matching  the  control  surface  deflection  and  rate  of  * 

*  deflection  each  time  a  new  command  is  issued  (initial  conditions  are  * 

*  all  taken  as  zero) .  Data  is  output  to  file  PL0T.DAT  every  * 

*  approximately  0.05  seconds.  * 


IMPLICIT  DOUBLE  PRECISION  (A-Z) 

INTEGER  I, J,MAXC0M,NC0M 
PARAMETER (  MAXCOM  =  20, 

+  PI  =  3. 1415926535897932385D0, 

+  RTOD  =  180. DO/PI, 

+  DT0R  =  1 . /RTOD  ) 

*  A  maximum  of  20  commands  in  a  series . 

DIMENSION  TC (MAXCOM) , DC (MAXCOM) 

*  TC:  the  time  at  which  a  command  is  issued. 

*  DC:  the  commanded  deflection  angle. 

*  Computer  dependent  numbers . 

COMMON  /CMPTR/  SMLNUM,SN100,T0L,TL1000,TL1L0G 
SMLNUM  =  l.E-5 
22  SMLNUM  =  0.5*SMLNUM 
A  =  SMLNUM+1 . 0 
IF  (A  .GT.  1.0)  GOTO  22 
SN100  =  SMLNUM* 100 . 

C  The  error  tolerance  level  (normally  set  by  the  calling  routine) . 
T0L  =  IE-5 

TL1000  =  MAX(T0L/1000 . ,SN100) 

TL1L0G  =  L0G(TL1000) 

DATA  DELC , AMP , DELTA/3*0 . 0/ 

*  Zero  initial  state  and  command  state  for  control  system. 
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*  Open  and  read  input  file  SIMU.INP. 

OPEN (50 , FILE= ’ SIMU . INP ’ , STATUS= ’ OLD ’ ) 

READ (50,*)  TEND 
READ ( 50 , * )  WMAX , DDMX , ZETA 
READ (50,*)  NCOM 
IF  (NCOM  .GT.  MAXCOM)  THEN 
WRITE(ITERM,4010)  MAXCOM, NCOM 
WRITE ( IOUTL ,4010)  MAXCOM, NCOM 

4010  FORMAT(/’  FATAL  ERROR:  Array  dimensioning  restricts  you  to  ’, 
+13,’  commands.’ /’  You  have  specified  ’,13,’  commands.’/) 

STOP 

ELSE  IF  (NCOM  .LE.  0.0)  THEN 
WRITE (ITERM, 4020)  NCOM 
WRITE (IOUTL, 4020)  NCOM 

4020  FORMAT (/’  FATAL  ERROR:  It  is  meaningless  to  specify  ’,13, 

+’  commands.’/) 

STOP 
END  IF 

READ (50,*)  (TC(I) ,DC(I) , 1=1 , NCOM) 

DO  100  1=1, NCOM 
100  DC (I)  =  DC(I) *DTOR 

IF  (TEND  .LT.  TC(1))  THEN 
WRITE (ITERM, 4030)  TC(1),TEND 
WRITE (IOUTL, 4030)  TC(1),TEND 

4030  FORMAT (/’  FATAL  ERROR:  The  integration  BEGINS  when  the  first  ’, 
+’ command  is  issued’/’  at  t  =’,f7.1,’,  but  you  have  specified  ’, 
+ ’ that  it  ENDS  at  t  =’ ,f 7. 1 , ’ . ’/) 

STOP 
END  IF 


T  =  TC(1) 
TSTART  =  T 


*  Data  is  written  to  specially  formatted  output  file  PLOT.DAT. 

OPEN (5 1 , FILE= ’ PLOT . DAT ’ , STATUS= ’ UNKNOWN ’ ) 

WRITE(51 , *)  1,0, 0,0,0 
WRITE (51 , *) 

WRITE (51 , *) 

WRITE (51 , *) 

WRITE (51 , *) 

WRITE (51 , *) 

1=1 

200  CONTINUE 

CALL  CNSYS2 (TC ( I ) , DC ( I ) , DELC , TO , AMP , BETA , OMEGA , WMAX , DDMX , ZETA) 

I  =  1+1 

IF  (I  .LE.  NCOM)  THEN 
TSTOP  =  min(TC(I) ,TEND) 
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ELSE 

TSTOP  =  TEND 
END  IF 

J  =  NINT ( (TSTOP-T) / 0.05) 

IF  (J  .LE.  5)  J=J+1 
DT  =  (TSTOP-T) /J 

300  CONTINUE 

CALL  CSDEFL (T , DELTA , DELC , TO , AMP , BETA , OMEGA , ZETA) 

IF  (T  .EQ.  TSTART)  THEN 
WRITE (51 , *)  0,T,DELTA*RT0D 
ELSE 

WRITE (51 , *)  1 , T,DELTA*RTOD 
END  IF 

T  =  T+DT 

IF  (T  .GT.  TEND+TOL)  GOTO  10 
IF  (T  .GT.  TSTOP+TOL)  GOTO  200 
GOTO  300 

10  CONTINUE 

WRITE(51 , *)  -1 
END 

SUBROUTINE  CNSYS2(T,D2,  DELC , TO , AMP , BETA , OMEGA ,  WMAX , DELDMX , ZETA) 


*  * 

*  (c)  Copyright  Her  Majesty  the  Queen  in  Right  of  Canada,  as  * 

*  represented  by  the  Department  of  National  Defence.  * 

*  * 


*  CNSYS2  accepts  a  new  commanded  deflection,  D2,  for  a  control  surface  * 

*  at  time  T  and  matches  a  second  order  response  to  the  current  control  * 

*  surface  position  and  speed  characteristics;  this  response  depends  on  * 


*  the  control  system  characteristics  WMAX,  DELDMX  and  ZETA.  The  * 

*  parameters  describing  the  response  to  the  previous  command,  DELC  * 

*  through  OMEGA,  are  upgraded  to  the  new  response  parameters.  * 

*  * 

*  T  -  Current  time.  * 

*  D2  -  Newly  commanded  deflection  angle,  effective  at  time  T.  * 

*  * 

*  DELC  -  Old  commanded  deflection  angle,  changed  to  D2  on  output.  * 

*  TO  -  Time  of  issuance  of  old  command,  changed  to  T  on  output.  * 

*  AMP  -  Oscillation  amplitude  of  response  at  time  previous  command  * 

*  was  issued  (T=T0) ,  output  value  is  for  new  command.  * 

*  BETA  -  Phase  shift  of  oscillation  of  previous  response,  output  * 

*  value  is  for  new  command.  * 

*  OMEGA  -  Frequency  at  which  control  system  reacts.  Input  value  is  * 

*  for  previous  command;  output  value  is  for  new  command.  * 

*  * 

*  Fixed  parameters:  * 
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*  WMAX  -  The  maxmimum  frequency,  omega,  at  which  the  control  system  * 

*  can  respond;  OMEGA  =  WMAX  unless  this  results  in  a  maximum  * 

*  rate  of  deflection  exceeding  DELDMX.  * 

*  DELDMX  -  The  maximum  rate  at  which  the  control  surface  is  capable  of  * 

*  deflecting;  a  positive  number.  * 

*  ZETA  -  The  damping  of  the  control  system.  * 

*  * 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

PARAMETER (  PI  =  3 . 1415926535897932D0 , 

+  RTOD  =  180. DO/PI, 

+  ZERO  =  0.0, 

+  ONE  =  1.0  ) 

C  Relative  error  used  in  this  subroutine  (TL1000)  is  1/1000  that  used  in 
C  the  calling  routine  so,  for  practical  purposes,  this  routine  is  exact. 
COMMON  /CMPTR/  SMLNUM,SN100,T0L,TL1000,TL1L0G 

SAVE  ZSAVE1 , ZSAVE2 , Z2 , RTZ , ACZ , ACZ2 , ZETRTZ , 

+  C0F7 , TEST , C0F6 , C0F5 , C0F4 , C0F3 , C0F2 

DATA  ZSAVE1 , ZSAVE2/2*-l . / 

C  Some  calculations  need  not  be  repeated  if  ZETA  doesn’t  change  from 
C  call  to  call. 

IF  (ZSAVE1  .NE.  ZETA)  THEN 
ZSAVE1  =  ZETA 

Z2  =  ZETA*ZETA 
RTZ  =  SQRT(1-Z2) 

ACZ  =  ACOS (ZETA) 

ACZ2  =  ACZ+ACZ 
ZETRTZ  =  ZETA/RTZ 

C0F7  =  -(Z2*(Z2*(Z2*7. 1168-26. 0928) -46 . 0728) -3 . 0375) /4354 . 56 
TEST  =  (TL1000/ C0F7) **0 . 142857 143 
END  IF 

50  CONTINUE 

IF  (AMP  .NE.  ZERO)  THEN 

WT  =  OMEGA* (T-TO) 

A1  =  -ZETA*WT 

C  If  original  amplitude  has  decayed  by  a  factor  of  TL1000,  AMP  =  0.0. 

IF  (A1  .GT.  TL1L0G)  THEN 
A1  =  AMP*EXP (Al) 

A2  =  RTZ*WT+BETA 
DELO  =  DELC-A1*SIN (A2) 

DELDO  =  0MEGA*A1*SIN (A2-ACZ) 

ELSE 

AMP  =  0.0 
GOTO  50 
END  IF 

C  If  DELDO  is  effectively  zero,  take  a  shortcut. 
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IF  (ABS(DELDO)  .LT.  TLIOOO)  THEN 
AMP  =  0.0 
DELC  =  DELO 
GOTO  50 
END  IF 

DELDEL  =  D2-DEL0 


C  Calculate  frequency  limited  characteristics. 

C  Evaluate  BETA  keeping  in  mind  that  DELDEL  may  be  zero . 

A1  =  ZETA*DELDEL-DELDO/WMAX 
BETA  =  ATAN2 (RTZ*ABS (DELDEL) , A1*SIGN (ONE , DELDEL) ) 

IF  (BETA  .GT.  0.7  .AND.  BETA  .LT.  2.4)  THEN 
AMP  =  DELDEL/SIN (BETA) 

ELSE 

AMP  =  Al/RTZ/ COS (BETA) 

END  IF 

DDMAX  =  AMP*WMAX*RTZ*EXP (ZETRTZ* (BETA-ACZ2) ) 

C  If  BETA  .LT.  ACZ2  and  if  the  maximum  rate  of  deflection  is  greater 
C  than  the  limit,  the  solution  is  rate  limited. 

IF  (ABS (DDMAX)  .GT.  DELDMX  .AND.  BETA  .LT.  ACZ2)  THEN 

DDMAX  =  SIGN (DELDMX, DELDEL) 

RATIO  =  DELDO /DDMAX 

C  See  if  Beta  can  be  solved  for  with  Taylor  series  or  if  Newton- 
C  Raphson  method  should  be  used. 

A1  =  MAX (ZERO, 1. -RATIO) 

EPS  =  SQRT (Al+Al) 

IF  (EPS  .LE.  TEST)  THEN 

C  Calculate  BETA  directly  using  its  Taylor  series  about  BETA  =  ACZ2. 
IF  (ZSAVE2  .NE.  ZETA)  THEN 
ZSAVE2  =  ZETA 

C0F6  =  ZETA*  (Z2*(Z2*1. 6-15. 6)-8.  D/850.5 
C0F5  =  (Z2* (Z2*6 . 4+62 .4) +8. 1)/1728. 

C0F4  =  -ZETA* (Z2*l . 6+2 . 7) /54 . 

C0F3  =  (Z2*8.+3.)/72. 

C0F2  =  -ZETA/3. 

END  IF 

BETA  =  ACZ2-RTZ*EPS* (EPS* (EPS* (EPS* (EPS* (EPS* 

+  (EPS*C0F7+C0F6) +C0F5) +C0F4) +C0F3) +C0F2) +1 . ) 

ELSE 

C  Calculate  BETA  using  the  Newton-Raphson  iteration.  Guess  an  initial 
C  value  for  BETA:  choose  initial  value  so  subsequent  iterations  march 
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C  BETA  monotonically  towards  its  correct  value  without  overshooting. 

IF  (RATIO  .GE.  0.0)  THEN 
C  a  minimum  value ! 

BETA  =  ACZ 
ELSE 

C  the  inflection  point  or  zero! 

BETA  =  MAX(ACZ2+ACZ-PI , ZERO) 

END  IF 

C  Enter  iterative  loop. 

A1  =  RATI0*RTZ 
100  A2  =  BETA-ACZ 

A3  =  A2-ACZ 

DB  =  RTZ/SIN (A3) * (SIN (A2) -A1*EXP (ZETRTZ*A3) ) 

BETA  =  BETA+DB 

C  Test  to  see  if  more  accuracy  is  required  (BETA  in  radians) . 

IF  (ABS(DB)  .GT.  TL1000)  GOTO  100 

END  IF 

C  The  rate  limited  BETA  is  now  known.  Calculate  the  rate  limited  OMEGA, 
C  saving  AMP  for  later  calculation. 

AMP  =  DELDEL/SIN (BETA) 

OMEGA  =  DDMAX/AMP/RTZ*EXP (ZETRTZ* (ACZ2-BETA) ) 

C  Note  that  DELDEL  cannot  be  zero;  if  it  is  the  solution  is 
C  necessarily  frequency  limited. 

ELSE 

C  The  frequency  limited  OMEGA. 

OMEGA  =  WMAX 

END  IF 

ELSE 

C  AMP  equal  zero,  a  special  case. 

DELDEL  =  D2-DELC 

IF  (ABS (DELDEL)  .LT.  TL1000)  THEN 
DELC  =  D2 
RETURN 
END  IF 

AMP  =  DELDEL/RTZ 
BETA  =  ACZ 

C  The  rate  limit  prediction  of  OMEGA: 

DDMAX  =  S I GN(DELDMX, DELDEL) 

OMEGA  =  DDMAX/DELDEL*EXP (ZETRTZ*ACZ) 

C  Convert  to  frequency  limit  prediction  if  appropriate. 

OMEGA  =  MIN (OMEGA, WMAX) 
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END  IF 


C  AMP  =  (D2-DELC) /SIN (BETA) ,  BETA,  and  OMEGA  have  now  been  determined. 

DELC  =  D2 
TO  =  T 

write(*,*)  ’  OMEGA, BETA  =  ’ , omega, beta*rtod 
write(*,*)  ’  DDMAX  =  ’ ,deldel/sin(beta) *omega*rtz 

+  *exp(zetrtz* (beta-acz2) ) 

END 

SUBROUTINE  CSDEFL (T , DELTA ,  DELC , TO , AMP , BETA , OMEGA , ZETA) 

*  Calculate  the  deflection,  DELTA,  of  a  control  surface  at  time  T.  The  * 


*  charateristics  of  the  second  order  control  system  for  the  control  * 

*  surface  are  given  by:  * 

*  * 

*  DELC  -  Commanded  deflection  angle.  * 

*  TO  -  Time  command  was  issued.  * 

*  AMP  -  Amplitude  of  oscillation  of  control  system  at  time  command  * 

*  was  issued.  * 

*  BETA  -  Phase  shift  of  oscillation.  * 

*  OMEGA  -  Response  frequency  at  which  oscillation  occurs.  * 

*  ZETA  -  Damping  of  control  system.  * 


IMPLICIT  DOUBLE  PRECISION  (A-Z) 

SAVE  ZSAVE , RTZ 

C  Relative  error  used  in  this  subroutine  (TL1000)  is  1/1000  that  used  in 
C  the  calling  routine  so,  for  practical  purposes,  this  routine  is  exact. 
COMMON  /CMPTR/  SMLNUM,SN100,T0L,TL1000,TL1L0G 
DATA  ZSAVE/- 1./ 

IF  (AMP  .EQ.  0.0)  THEN 
DELTA  =  DELC 
ELSE 

WT  =  OMEGA* (T-TO) 

ZWT  =  -ZETA*WT 

C  If  original  amplitude  has  decayed  by  a  factor  of  TL1000,  DELTA  =  DELC. 
IF  (ZWT  .GT.  TL1L0G)  THEN 
IF  (ZETA  .NE.  ZSAVE)  THEN 
RTZ  =  SQRT(1 . -ZETA*ZETA) 

ZSAVE  =  ZETA 
END  IF 

DELTA  =  DELC-AMP*EXP (ZWT) *SIN (RTZ*WT+BETA) 

ELSE 

DELTA  =  DELC 
END  IF 
END  IF 

END 
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Appendix  C:  Simulation  Plots 

Each  of  the  nine  pages  following  this  one  contain  graphical  output  from  a  rising  simulation 
using  the  coefficient  based  hydrodynamic  model.  The  first  simulation,  Figure  Cl,  forms  the 
basis  for  Simulations  2  through  8  (Figures  C2  through  C8).  Simulation  9  (Figure  C9)  is  based 
on  Simulation  8.  Hereafter,  Simulation  n  is  referred  to  as  Sn. 

To  simplify  comparisons,  the  scales  on  graphs  of  the  same  type  are  the  same.  The  only 
exceptions  are  the  (b)  and/or  (d)  plots  in  S6,  S7,  and  S9;  however,  these  contain  horizontal 
dashed  lines  showing  the  normal  location  of  the  limits. 

All  variables  plotted  are  listed  along  the  ordinate  axis  in  an  order  that  determines  the  color 
of  the  curve  showing  that  quantity  in  the  plot: 


Variable  1 

Red 

2 

Green 

3 

Blue 

4 

Black 

5 

Cyan  (light  blue) 

6 

Orange 

The  exceptions  are  the  Vai/VTi  curves  in  parts  (e)  of  the  figures.  These  follow  the  same  color 
scheme  used  in  Figure  8  in  the  main  text: 

val/vT1 

Red 

Va2/VT2 

Green 

Va3/VT3 

Cyan 

Va4/VT  4 

Blue 

and  are  easily  identified  as  the  only  curves 

in  parts  (e)  which  attain  a  value  of  1.0. 

All  the  plots  with  time  along  the  abscissa  axis  have  vertical  black  lines  marking  significant 

events.  Fine  dashed  vertical  lines  indicate  times  at  which  new  sternplane  commands  are  is- 

sued.  Fine  solid  vertical  lines  mark  the  times  when  the  first  and  fourth  MBTs  suddenly  empty 
(Vai/VTl  =  1);  the  fourth  tank  often  does  not  empty  before  the  boat  reaches  the  surface.  These 
events  are  also  marked  in  parts  (c)  of  the  figures  (the  trajectory)  by  the  plus  (+)  symbols. 

In  addition,  a  dark  vertical  line  with  circles  (o)  at  its  end  points  marks  the  time  at  which 
the  boat  becomes  unstable  in  roll.  This  is  determined  by  when  Us  =  0  in  parts  (a).  In  parts  (c), 
this  point  is  marked  by  just  the  circle  symbol. 

The  plots  end  when  the  boat  emerges  through  the  surface,  at  t  =  te.  This  is  defined  as 
the  time  at  which  the  top  of  the  forward  most  MBT  is  at  the  ocean  surface.  This  means  the 
simulations  finish  with  the  body  fixed  axes  origin  at  various  depths  (cf,  parts(c)  of  S6  and  S7) 
because  of  different  surfacing  pitch  angles. 

Parts  (f)  of  the  figures  show  0  and  4>.  As  explained  in  the  main  text,  can  jump  by 
7 r  when  w  changes  sign  (depending  on  v ),  which  it  does  in  the  first  10  seconds  or  so  of  SI 
and  comes  close  to  doing  in  several  of  the  other  simulations.  When  this  happens,  something 
approaching  a  discontinuity  appears  in  0;  it  would  be  a  discontinuity  if  v  were  exactly  zero 
when  w  changed  sign. 
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a)  Velocity  and  propulsion  states. 
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b)  Rotational  velocity  states. 


e) 


Figure  Cl  Simulation  1:  The  baseline  simulation.  Start  from  straight  and  level  flight  at  100  m 
depth  and  a  speed  of  3  m/s.  At  t  =  0,  command  a  sternplane  deflection  of  -20  degrees  to  pitch 
up.  Simultaneously  blow  ballast  normally  and  command  a  speed  increase  to  6  m/s.  Adjust  the 
sternplane  deflection  to  limit  the  pitch  angle  to  20  degrees. 
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a)  Velocity  and  propulsion  states. 
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b)  Rotational  velocity  states. 


t  (s) 

d)  Euler  angles. 


e)  Mass  and  blow  characteristics. 


0  10  20  30  40  50  60  70 


t  (s) 

f)  Incidence  and  sternplane  deflection. 


Figure  C2  Simulation  2:  The  same  as  SI  except  commanded  speed  is  not  increased. 
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e)  Mass  and  blow  characteristics. 


f)  Incidence  and  sternplane  deflection. 


Figure  C3  Simulation  3:  The  same  as  SI  except  the  initial  speed  is  4.5  m/s,  increased  to 
9  m/s  beginning  at  t  =  0;  that  is,  the  speeds  are  increased  by  50%,  approximately  doubling  the 
dynamic  pressure. 
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a)  Velocity  and  propulsion  states. 
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b)  Rotational  velocity  states. 
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d)  Euler  angles. 


Figure  C4  Simulation  4:  The  same  as  SI  except  sternplane  control  is  not  used  to  pitch  the 
nose  up. 
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a)  Velocity  and  propulsion  states. 
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b)  Rotational  velocity  states. 


t  (s) 

e)  Mass  and  blow  characteristics. 
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f)  Incidence  and  sternplane  deflection. 


Figure  C5  Simulation  5:  The  same  as  SI  except  the  MBTs  are  blown  using  the  emergency  air 
supply  which  provides  twice  the  mass  of  air  as  for  normal  blowing. 
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a)  Velocity  and  propulsion  states.  b)  Rotational  velocity  states. 


Figure  C6  Simulation  6:  The  same  as  SI  except  the  sternplanes  are  used  to  curtail  the  pitch 
(once  50  m  depth  is  achieved)  so  the  boat  surfaces  with  a  pitch  angle  under  10  degrees. 
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a)  Velocity  and  propulsion  states. 


b)  Rotational  velocity  states. 


e)  Mass  and  blow  characteristics. 


f)  Incidence  and  sternplane  deflection. 


Figure  C7  Simulation  7:  The  same  as  SI  except  the  pitch  angle  is  allowed  to  increase  to  35 
degrees. 
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a)  Velocity  and  propulsion  states. 
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b)  Rotational  velocity  states. 


e)  Mass  and  blow  characteristics. 


f)  Incidence  and  sternplane  deflection. 


Figure  C8  Simulation  8:  The  same  as  SI  except  the  initial  roll  angle  is  increased  to  2  degrees 
by  setting  yGo  (normally  0)  to  -11  cm. 
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a)  Velocity  and  propulsion  states. 
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f)  Incidence  and  sternplane  deflection. 


Figure  C9  Simulation  9:  The  same  as  S8  except  the  sternplanes  are  used  to  limit  the  pitch 
angle  to  10  degrees  throughout  the  maneuver. 
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